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Chapter 1 
Introduction 



The central objective of particle physics is to study the basic building blocks of matter, and 
to explore the way they bind together and interact with each other. Based on the interaction 
type one can distinguish between four different forces acting on the elementary particles: these 
are the gravitational, electromagnetic, weak and strong interactions. While a proper quantum 
description of gravity is not yet established, the latter three interactions are summarized in 
the structure which is called the Standard Model. The sector of the Standard Model that 
covers the strong force is described by a theory called Quantum Chromodynamics (QCD). The 
elementary particles of QCD notably differ from those of e.g. the electromagnetic interaction: 
they cannot be observed directly in nature. These particles - quarks and gluons - only show 
up as constituents of hadrons like the proton or the neutron. 

On the other hand, according to QCD, in certain situations quarks are no longer confined 
inside hadrons. One of the most important properties of QCD is asymptotic freedom, which 
implies that the interaction between quarks vanishes at asymptotically high energies. Due 
to asymptotic freedom, under extreme circumstances - namely, at very high temperature or 
density - quarks can be liberated from confinement. At this point a plasma of quarks and 
gluons comes to life that is substantially different as compared to the system of confined 
hadrons. This plasma state is referred to as the quark-gluon plasma (QGP). 

1.1 The quark-gluon plasma 

There are two situations in which QGP is believed to exist (or have existed): first in the early 
Universe and second in heavy ion collisions. In both cases the temperature - and with it the 
average kinetic energy of quarks - is high enough to overcome the confining potential that is 
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present inside hadrons and due to this the very different plasma state can be formed. For the 
case of the early Universe this plasma state was realized until about 10~^ seconds after the Big 
Bang, when the temperature sank below a critical temperature Tc ~ 200 MeV^. The properties 
of the transition between the 'hot', deconfined QGP and the 'cold', confined hadrons play a 
very important role in the understanding of the early Universe [1]. The two above forms of 
strongly interacting matter can be thought of as phases in the sense that the dominant degrees 
of freedom in them are very different: colorless^ hadrons in the cold, and colored objects in 
the hot phase. In accordance with this, the transition can be treated as a phase transition in 
the statistical physical sense. 

One of the most important properties of the transition is its nature. We define a transition 
to be of first order if there is a discontinuity in the first derivative of the thermodynamic 
potential. For a second- order transition there is a jump in the second derivative, i.e. the first 
derivative is not continuously differentiable. For an analytic transition no such singularities 
occur, one may refer to such a process as being a crossover. In the case of a first-order phase 
transition - like the boiling of water - the quark-gluon plasma 
would reach a super-cooled state in which smaller bubbles of the 
favored, cold phase can appear. As the system aspires towards 
the minimum of the free energy, the large enough (supercritical) 
bubbles can further grow and after a while merge with each 
other (see illustration in figure 1.1, taken from [2]). This process 
can be treated as a jump through a potential barrier from a local 
minimum to a deeper minimum; from a so-called false vacuum 
to the real vacuum. On the other hand the transition can also 
be continuous (second order or crossover) - in this case no such 
bubbles are created and the transition between the two phases 
occurs uniformly. 

The cosmological significance of the above phenomenon is 
that if such bubbles indeed appeared then at the phase bound- 
aries specific reactions can take place that one would be able 
to observe in the cosmic radiation. Such a transition may also 
have a strong effect on nucleosynthesis. The boundaries of the 
bubbles can furthermore collide and as a result produce gravitational waves that may also be 

^This is equivalent to about 10^^ degrees Celsius. 

^The analogue of charge in QCD is called color, see section 2.1. 
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detected [3]. According to lattice calculations of QCD however, the transition from QGP to 
confined, hadronic matter is much likely to be an analytic crossover [4]. This means that the 
transition goes down continuously as illustrated in figure 1.2. 



Hot 



Aging of the Universe 



Cold 



Figure 1.2: Illustration of a crossover transition. Mesons (quark-antiquark bounded states) and 
hadrons (three-quark bounded states) appear continuously while the Universe cools down to its cold 
phase. 



The QCD transition plays a very relevant role also from the point of view of heavy ion 
collisions. It is now widely accepted that in high energy collisions of heavy ions conditions re- 
sembling the early Universe can be generated and the plasma phase of quarks can be recreated. 
Recently, in a collision of gold nuclei at the Relativistic Heavy Ion Collider (RHIC), an initial 
temperature beyond 200 MeV was reached [5]. There are also further signals indicating that 
the QGP has indeed been created. One of these signals is the phenomenon of jet quenching. A 
jet is a beam of secondary particles that originates from the high-momentum quark that was 
broken out of the incoming protons or neutrons. Interactions between the jet and the hot, dense 
medium produced in the collision are expected to lead to a loss of the jet energy. Evidence for 
jet quenching has indeed been found at the Relativistic Heavy Ion Collider (RHIC) [6]. 



1.2 The phase diagram of QCD 

Just like at high temperature T, also at large quark densities p we expect (in agreement with 
asjTiiptotic freedom) the coupling between quarks to decrease and the QGP phase to be created. 
In the statistical physical approach to QCD thermodynamics, the net density of quarks^ can be 

^The density of quarks minus that of antiquarks. 



3 



1. INTRODUCTION 



controlled by a chemical potential Zero chemical potential corresponds to a situation where 
the density of quarks and antiquarks is the same. The phases of strongly interacting matter 
can be represented ona/i — Torp — T phase diagram (see a possible depiction on figure 1.3, 
taken from [7]). On this diagram phases are separated by transition lines that can represent 
either first-order or second-order phase transitions or continuous transitions (crossovers). The 
T, /i parameters during the cooling down of the early Universe or a high energy collision also 
draw a trajectory on the phase diagram. These trajectories are contained in the small /i region 
of the phase diagram, since in both cases the number of quarks and antiquarks are roughly 
the same. Accordingly, this situation represents a thermodynamic system having zero or small 
chemical potential. 

By increasing the density and keeping the temperature fixed - i.e. compressing hadronic 
matter - one can also move into the deconfined phase of quarks. We know far less about this 
region of the phase diagram, which is thought to exhibit phenomena like color flavor locking 
or color superconductivity. On the other hand, the low fi area is much better understood 
and theoretically more tractable. Besides phenomenological interest, the detailed structure of 
this area (like the transition temperature or the curvature of the transition line) is relevant 
for contemporary and upcoming heavy ion experiments. While most of the ongoing experi- 
ments like those at LHC or RHIC concentrate on achieving very high energies and thus small 
chemical potentials, there are projects that aim for regions of the phase diagram with larger 
densities (RHIC II, FAIR)^. Designing these next generation experiments can benefit greatly 
from developing theoretical understanding of the phase diagram. 

According to lattice simulations, at zero chemical potential the transition is an analytic 
crossover [4]. This is represented in figure 1.3 by the smooth transition between the white 
and the orange regions. A possible scenario about the p > region of the phase diagram 
is that a first-order transition emerges (yellow band), which also implies the existence of a 
critical endpoint (orange dot). Such a critical endpoint corresponds to a second-order phase 
transition that belongs to a given universality class. It could also happen that the transition 
is continuous also at larger fi values, and then no critical endpoint exists. There have been 
lattice indications favoring both scenarios, so this is still an open question. See e.g. [8] arguing 
for, and [9] against the existence of a critical endpoint. 

A further important characteristic of the system is its equation of state as a function of the 
temperature. The equation of state (EoS) is also sensitive to the transition between hadronic 

^In a heavy ion collision the density of the system is controlled by the center of mass energy y^s^vtv and 
the centrality of the collision. 
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Net Baryon Density 



Figure 1.3: A possible depiction of the phase diagram of QCD in the space of the state parameters: 
the temperature T and the baryon density (which equals three times the quark density). Phases are 
separated by a crossover transition at zero chemical potential. At larger values of p a first-order phase 
transition may take place, which is indicated by the yellow band. The crossover and first-order lines 
must then be separated by a critical endpoint. Regions at large densities and small temperatures are 
thought to describe the interior of dense neutron stars. At even higher densities exotic phases like a 
color superconductor are expected. Figure taken from [7]. 

phase and the QGP and thus also plays a very important role in high energy particle physics. 
Moreover, recent results from RHIC imply that the high temperature quark-gluon plasma 
exhibits collective flow phenomena. It is also conjectured that the description of hot matter 
under these extreme circumstances can be given by relativistic hydrodynamic models. In turn, 
these models depend rather strongly on the relationship between thermodynamic observables, 
summarized by the equation of state. The EoS can be calculated using perturbative methods, 
but unfortunately, such expansions usually converge only at temperatures much higher than 
the transition temperature. Therefore the lattice approach (as a non-perturbative method) is 
a suitable candidate to study the EoS in the transition region T ~ T^. 

1.3 Structure and overview 

In this thesis I concentrate on the low /i, high T region of the phase diagram, which - according 
to the above remarks - is interesting for both the context of the evolution of the early Universe 
and heavy ion collisions. 

The thesis is structured as follows. First I present a brief introduction to the theoretical 
study of the QGP. This includes the definition of the underlying theory, QCD, and the method 
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with which QCD can be represented and studied using a finite discretization of the variables on 
a four-dimensional lattice (see chapter 2). Using the lattice formulation one can study various 
thermodynamic properties of strongly interacting matter. This is investigated in detail for the 
case of vanishing quark density and also for the case where a positive net quark number is 
present. The analysis of the latter system entails a fundamental problem, so separate sections 
are devoted to this issue (see chapter 3). 

After having discussed some of the basic elements of lattice QCD thermodynamics, I turn to 
present several new results regarding the transition between confined hadrons and the quark- 
gluon plasma. These results are divided into three separate chapters. First, in chapter 4 
I present the study of the phase diagram at small chemical potentials. In this project the 
pseudocritical temperature and the nature of the QCD transition are analyzed as a function 
of the quark density. With the help of these functions the phase diagram of QCD can be 
reconstructed. In particular, the curvature of the transition line lying between the two phases 
is determined, and the possibility of the existence of a critical endpoint is also addressed. 
Preliminary results regarding the curvature have been published in [10], while the full result 
has been published recently [11]. This work was done in collaboration with Zoltan Fodor, 
Sandor Katz and Kalman Szabo. My contributions to the project were the following: 

• I have developed and implemented a method to define the curvature without the need to 
fit the fi > data. This definition also gives information regarding the relative change 
in the strength of the transition. 

• I have performed all of the simulations and measured the Taylor-coefficients necessary 
for this definition. By means of a multifit to data measured at various lattice spacings 
I determined the curvature of the transition line separating the confined and deconfined 
phase. 

Afterwards I turn to show results regarding the equation of state. The central quantity 
here is the pressure as a function of the temperature. After a brief overview of the litera- 
ture and a discussion about how one can determine the pressure on the lattice I present the 
results. In chapter 5 I study the EoS with dynamical quarks; this work has been recently pub- 
lished [12]. I participated in this project as member of the Budapest- Wuppertal collaboration. 
My contributions were: 

• I have developed a multidimensional integration scheme that can be used to reconstruct a 
smooth hypersurface using scattered gradient data. I used this approach to determine the 
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pressure of QCD in tlie two-dimensional parameter space spanned by the gauge coupling 
and the light quark mass. In this approach it is straightforward to study the quark mass 
dependence of the EoS and to estimate the systematic error in the pressure. 

• I have measured the charm condensate on part of the dynamical configurations and eval- 
uated the charm contribution to the pressure and to other thermodynamic observables. 

The multidimensional integration method is applicable on a general level and thus was also 
published individually [13]. The method is summarized in more detail in appendix A. 

Finally, in chapter 6 I show results regarding the high temperature EoS in the pure glu- 
onic theory. At these previously unreached temperatures it becomes possible to carry out a 
comparison to improved or resummed perturbation theory. Furthermore, the non-perturbative 
contribution to the trace anomaly is also quantified. This is also a project of the Budapest- 
Wuppertal collaboration, currently under publication [14], with preliminary results already 
published earlier [15]. My contributions were the following: 

• I have evaluated the pressure using a multi-spline fit to data measured at various lattice 
spacings and extracted the continuum extrapolated curve from this fit. 

• I have compared and matched lattice results with improved and resummed perturba- 
tive formulae. The free parameters were fitted to best reproduce the lattice results at 
high temperature. Furthermore, I quantified the non-perturbative contribution to the 
interaction measure with either a constant or a logarithmic ansatz. 
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Chapter 2 



Quantum Chromodynamics and the 
lattice approach 

2.1 A theory for the strong interactions 

Nowadays it is widely accepted that QCD is the appropriate tool to treat the interactions 
between quarks and gluons. For a long time it was however unclear what kind of theory could 
describe the strong forces [16]. The electron-hadron collision experiments of the late sixties 
- namely, the Bjorken-scaling of the structure functions in such scatterings - suggested that 
electrons scatter off almost free, point-like constituents. The Bjorken-scaling implies that these 
constituents - the quarks - should interact weaker at shorter distances (at larger energies). 
Since for electromagnetism an appropriate describing theory was Quantum Electrodynamics 
(QED), it was reasonable to search for another quantum field theory (QFT) that succeeds 
to describe the dynamics of quarks. While most QFTs fail to fulfill the above requirement, 
it was proven in 1973 that non-Abelian gauge field theories on the other hand possess this 
property of exhibiting a weaker force at shorter distances. This property was named asymptotic 
freedom [17, 18]. 

QCD is a non-Abelian quantum field theory, which differs from its Abelian relative QED 
in the fact that here the symmetry transformations of the underlying theory cannot be inter- 
changed. In other words, the corresponding symmetry consists of non-commutative generators: 
the gauge group here is SU(3) instead of U(l). Not much after the study of Gross, Wilczek 
and Politzer it was proven that non-Abelian gauge theories are not just a possible candidate 
to play the role of the theory of strong interactions; they represent the only class of theories 
in four space-time dimensions that exhibit asymptotic freedom. 
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The new symmetry (the gauge symmetry) that is generated by the non-commutative algebra 
describes a new type of property that quarks - compared to electrons - possess. This property 
was named color due to the apparent analogy of the system to color mixing: while quark 
fields carry a color quantum number, three of them can only build up a proton if the resulting 
combination possesses no color indices. In other words, three colored quarks can only be 
confined inside a proton if their 'mixture' is colorless. If one assumes that only such colorless 
states can be realized, then it is obvious that quarks as isolated particles cannot exist. This 
property of QCD is called confinement - referring to the fact that three quarks are confined 
inside a proton. 

Just like in QED, in QCD there is also a mediating particle corresponding to the gauge field 
that transmits the interaction between color-charged quarks. The analogon of the QED photon 
is the gluon. While the photon does not have an electric charge (due to the non- commutative 
nature of the gauge group), the gluon possesses a color quantum number. Because of this, 
contrary to the photon, the gluon also couples to itself; this is the reason for the fact that 
instead of the 1/r-like decaying potential of QED, in QCD a linearly rising potential appears 
between color-charged objects. The linear potential can be interpreted as a spring that binds 
quarks to each other inside a proton. This is just another way to describe the phenomenon of 
confinement. 

In accordance with the above remarks, the Lagrangian density of QCD contains the fermion 
and antifermion fields ip and ^/^ ("quarks" and "antiquarks" ) and the gauge field ("gluons"). 
The dynamics of the gauge field is governed by the field strength tensor 

= ^^,A, - d,A^ - ig[A^, A,] (2.1) 

Furthermore, the interaction between gluons and quarks is determined by the minimal coupling, 
as contained in the covariant derivative 

D^ = d^- igA^ (2.2) 

Putting all this together, the QCD Lagrangian describing Nq number of different quark flavors 
(in Minkowski space-time) is: 

1 

L = --Tr [F^^Fn + 5Z (^7m^^ - m,) (2.3) 

9=1 

The parameters of the theory are the gauge coupling g and the masses of the quarks niq. In 
nature there are six quark flavors (up, down, strange, charm, bottom and top) and thus six 
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masses. However, the contribution of heavy quarks is usually negligible at the non-perturbative 
scale of Tc ~ 200 MeV and only the first three or four quark species play a significant role. 
Furthermore, the difference between the up and down masses is very small (compared to T^) 
and thus it is a good approximation to take = = rriud- 

QCD is an SU{3) gauge theory, i.e. the Lagrangian (2.3) is invariant under SU{3) trans- 
formations. Quarks and antiquarks transform according to the fundamental representation of 
the gauge group, on the other hand gluons are placed in the adjoint representation. Accord- 
ingly, the fermionic fields have three color components and gluons contain eight degrees of 
freedom. Furthermore, in view of the Lorentz group quarks are bispinors and thus have four 
spin-components: 

quarks: ipq" c = 1 . . . 3, a = 1 ... 4, q = u,d, s,c,t,b (2.4) 

8 

gluons: A^ = Y,KT^ (2-5) 

a=l 

with Ta being the infinitesimal generators of the gauge group, usually represented by the eight 
Gell-Mann matrices. The Lorentz- and color indices of the quark field will be suppressed in the 
following. Moreover, throughout the thesis the different quark flavors will be identified using 
their first letter in the index of the field e.g. ipu will stand for an up quark. 

2.2 Perturbative and non-perturbative approaches 

According to the property of asymptotic freedom, for short-range reactions the interaction is 
weak and thus such processes can be safely analyzed by perturbation theory. On the other 
hand the temperature scale on which perturbative expansions converge is extremely high, and 
therefore in order to study e.g. the transition to QGP it is necessary to assess the dynamics of 
quarks using a non-perturbative approach. In the second half of the 20th century a new theory 
was constructed that is based on mathematical concepts and can be investigated through 
numerical simulations: lattice gauge theory. The main idea of this approach is to restrict the 
fields of the QCD Lagrangian to the points of a four-dimensional lattice. This theory is the only 
non-perturbative, systematically adjustable approach to QCD. Through lattice gauge theory 
we can gain information about quarks and gluons solely using first principles - the Lagrangian 
of QCD. 

The appearance of ultraviolet divergences, being a familiar property of quantum field theo- 
ries, is also present in QCD. A possible way to deal with such infinities in QFT is to introduce 
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some type of regularization (e.g. a cutoff) that makes tlie divergent Feynman-amplitudes 
mathematically tractable. After a proper renormalization of these amplitudes, the next step 
is the removal of the regularizing constraint. Necessarily, due to the redefinition of the renor- 
malized quantities - which will then converge to a finite value as the regularization is removed 
- the bare parameters of the theory will also become cutoff-dependent. The renormalization 
program as implemented on the lattice will be discussed in details in section 3.2. 

This procedure of regularization and renormalization can be realized using lattice gauge 
theories in an instructive manner. The lattice itself plays the role of the regulator, since the 
introduction of a finite lattice spacing a is equivalent to setting a cutoff 27r/a in momentum- 
space. This way on a lattice of finite size obviously every amplitude is going to be finite. The 
removal of the finite lattice spacing is called the continuum limit: this is usually done through 
some kind of extrapolation to a — )■ 0. 

The transition from confined hadrons to QGP is definitely a phenomenon that is only 
accessible by the lattice approach, since the relevant temperature scale of Tc ~ 200 MeV is still 
very far from the perturbative region. In the following sections I will address basic elements 
of lattice QCD, starting from the discretization of the variables, for the gauge fields and also 
for the fermion fields. The latter entails a fundamental problem that is stated in the so-called 
Nielsen- Ninomiya no-go theorem. Among the possible discretizations I will investigate the 
staggered version, since this was used to obtain all the results that are presented in this thesis. 

2.3 Quantum field theory on the lattice 

In this section a very brief overview of lattice gauge theory is given. A full and detailed 
introduction can be found in e.g. [19], [20] or [21]. 

As mentioned above, the lattice can be thought of as a possible regulator of the divergent 
Feynman-amplitudes. The lattice approach can be formulated through the functional-integral 
formalism, which is a generalization of the quantum mechanical path integral of Feynman. 
This formalism readily shows how it becomes possible to treat the system non-perturbatively, 
without any use of perturbation theory. 

According to the path integral in quantum mechanics, the Green's function of a particle 
propagating from qi to q2 in the time interval [ti, ^2] can be written as an integral over various 
possible paths: 




(2.6) 
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where S[q] is the action that belongs to the path given by q. This expression tells us to take 
every possible classical path that starts from qi and ends at q2, and sum the corresponding 
phase factors in the integrand. These factors, however, oscillate very strongly, so the calculation 
has to be carried out in Euclidean space-time instead of the usual Minkowski space-time, which 
can be achieved using the t — > —ir substitution. This latter change of variables is referred 
to as a Wick-rotation, after which time formally flows in the direction of the imaginary axis^. 
Now in the above expression the exponent of the integrand changes to minus the Euclidean 
action Se, which is the (imaginary) time integral of the Euclidean Lagrangian Le'- 



The above formula is mathematically well defined if one divides the time interval into pieces, 
calculates the finite sum for a given and then takes the N ^ oo limit. Expression (2.7) 
should be interpreted in this sense. 

The same procedure can be carried out in quantum field theories also. Here, instead of 
a finite number of degrees of freedom one deals with a field at each spacetime-point. While 
in quantum mechanics all information about the system is contained in the Green's func- 
tion (2.6), in a scalar QFT this role is played by an infinite set of ground state expectation 
values {^p{xi) . . . ip{xn))- In the same manner as before one obtains 



Here, analogously to the quantum mechanical case, the Tiip symbol indicates that the integrand 
has to be evaluated for each possible field configuration. Again, this expression is only well- 
defined for a finite number of degrees of freedom, so now not just time, but space also has to be 
discretized. This means that the field (f has to be restricted to the sites of a four-dimensional 
lattice. The result for the above functional-integral is given by taking the limit where the 
discretization is taken infinitely fine. 

More interesting from the quark-gluon point of view is the case of gauge fields and fermionic 
fields, since these appear in the Lagrangian of QCD shown in (2.3). Let us consider the case 

^Iii order to interpret results obtained after the Wick- rotation a continuation back to real time is of course 
due. However, for time-independent processes this is not necessary. 




(2.7) 




(2.8) 
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of only one flavor of quarks with mass m and coupling g. The latter usually enters the action 
in the combination /3 = G/g"^: 

Z = j DA^'DtpV^Jje-^^^^^'^'^''^''^^ (2.9) 

Since quarks are fermions, according to the spin-statistics theorem the fields ip and ip can be 
represented by anticommuting Grassmann-variables. With the functional integral given the 
next step is to put the theory on the lattice. However, the way the fields A^, ip and ip are 
discretized is not at all unique. Before I turn to the discretization definitions, it is instructive 
to interpret expression (2.9) a bit more thoroughly. 

2.3.1 Statistical physical interpretation 

If the action Se is bounded from below, the expression (2.9) - with its lattice discretized defi- 
nition - has the same form as the partition function of a classical statistical physical ensemble 
in four dimensions. This formal correspondence is valid at zero temperature. The quantum 
partition function Z = e~^^^ at a finite temperature T is on the other hand represented by 
a functional integral in which the integral of ZL^; in the imaginary time direction is restricted 
to a finite interval of length For bosonic fields periodic, for fermionic fields antiperiodic 

boundary conditions have to be prescribed in this direction. 

This interpretation justifies the lattice approach as a non-perturbative method as statistical 
physical methods can be used to calculate Green's functions like the one in (2.8). According to 
this analogy Z is called partition function and the Green's functions derived from Z are called 
correlation functions. The expectation value of an arbitrary operator can then be written as 

(0) = ly DA^2)^D^0e^^^[^'"^'^''"'^] (2.10) 

An important remark to make here is that while in the quantum theory the temperature is 
determined according to the Boltzmann-factors e~^^'^ , here T is proportional to the inverse of 
the size of the system in the Euclidean time-direction. In particular, on a lattice with lattice 
spacing a the temperature and the volume of the system are accordingly given as 

T = ^, V={Nsa)' (2.11) 

where Ng (iVj) is the number of lattice sites in the spatial (temporal) direction. In a usual 
computation one uses an A^^ x Nf lattice, so the sizes in the spatial directions are the same. It 

^The Boltzmann-constant is set here to unity: = 1- 
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is important to note that lattices where A^^ ^ Nt correspond according to (2.11) to a system 
with finite nonzero temperature; on the other hand lattices with A^^ > A^^ realize systems with 
roughly zero temperature. Also, the total volume of the system is given by V^d = V/T. 

2.3.2 Gauge fields and fermionic fields on the lattice 

As part of the regularization process we have to discretize the action Se, which is given by 
the four dimensional integral of the Euclidean Lagrangian density^ of QCD. It can be proven 
that in order to preserve local gauge invariance - which is of course indispensable to formulate 
a gauge theory - gauge fields must be introduced on the links connecting the sites rather 
then on the sites themselves (as in the case of the scalar field). Only this way are we able 
to represent local gauge transformations in a manner that fits the definition of the continuum 
transformation. Assigning the gauge fields of the QCD Lagrangian to the sites of the lattice 
would make this compliance impossible. On the links the gauge fields can be represented with 
Ufj, = 6*°^'^'' G SU{3) matrices^. This also implies that the Hermitian conjugate (i.e. the 
inverse) of the matrix representing a given link equals the matrix corresponding to the link 
pointing to the opposite direction: 



Here /t denotes the unit vector in the /i direction and n is the lattice site. Due to the trans- 
formation properties of ?7^, the simplest gauge invariant combination of gauge fields on the 
lattice can be constructed by taking the product of the links that build up a square (in the 
/i, u plane) 



and then calculating the trace of this expression. This trace is called the plaquette, based on 
which the action corresponding to pure gauge theory can be constructed. The resulting sum 
is the simplest real and gauge invariant expression that can be built using only gauge fields: 



The sum extends to every possible {fi — u) square on the four-dimensional lattice. The com- 
bination (2.14) is called the Wilson gauge action. Choosing f3 = 6/ it is straightforward 

^In the following the subscript E is dropped. 

^This combination ensures that as a — > the correct continuum theory is approached. 




(2.12) 




(2.13) 




(2.14) 
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to show that in the a — hmit the Wilson action approaches the continuum gauge action, 
namely the first term in (2.3).^ 

To obtain the total action of QCD one also has to take into account the fermionic contri- 
bution. According to the transformation rules of the fermionic fields ip and ip (which live on 
the sites of the lattice) other types of invariant combinations can also be composed. Since the 
QCD Lagrangian contains fermions quadratically, the general form of the action is written as 

S{U, ij, ij) = Sg{U) - ^M{U)ij (2.15) 

where M{U) is the fermion matrix. The elements of this 12A^ x 12A^ matrix (with being 
the number of lattice sites and 12 = 3 ■ 4 the number of colors times the number of Dirac- 
components) can be read off from the Lagrangian. The fermion matrix can be divided into 
a massless Dirac operator and a mass term: M = Jj) + ml. In the partition function the 
integration over the fermions can be analytically performed and gives the following result^: 

Z = y ^Df/DT/'DT/; e"^^^^)-'^*^^^)^ = j DU e~^^^^^ det M{U) (2.16) 

Here the integration measure for the fermions takes the simple product form 

D^ = Y[dtlin) (2.17) 

n 

the one over the gauge fields on the other hand depends on the 8 real parameters of the SU (3) 
group, and the integral has to be performed over the whole group. If the parameters on the 
^th link are denoted by a", the measure can be written as 

8 

2)?7 = JJ J(a^) JJ rfa^ (2.18) 

e a=l 

where the structure of the J{ag) Jacobi-matrix can be determined by requiring gauge invari- 
ance. This integration measure is called the Haar-measure. 

^Note that the Wilson action can be extended by an arbitrary term that vanishes as a — )• and the resulting 
action still converges to the same continuum action. I get back to this in section 2.4. 
^Recall that 4' and -ip are Grassmann- variables. 
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2.3.3 Fermionic actions 

The naive discretization of the fermionic part of the Lagrangian (2.3) fails to give the correct 
continuum limit even in the free case. Namely, the naive fermionic action 



^naive ^ ^ 



(2.19) 

gives a propagator for the free theory (where U = 1) which possesses not 1 but 16 poles in the 
lattice Brillouin-zone — vr/a < p < tt /a. Thus, the action (2.19) describes 16 quarks and does 
not converge to the continuum action as a — 0. This is a consequence of the Nielsen-Ninomiya 
theorem which states that this 'doubling' problem cannot be solved without breaking the chiral 
symmetry of the QCD action in the m — > limit. For the massless Dirac operator Ip, chiral 
symmetry means that 

{0,75} = 075 + 75^ = (2.20) 

is satisfied. In the continuum theory, although this symmetry would imply the conservation 
of an axial- vector current, the corresponding current has an anomalous divergence due to 
quantum fluctuations. On a lattice with finite lattice spacing however, this current is indeed 
conserved, and the corresponding extra excitations are just the above mentioned 'doublers'. 

The two most popular methods to circumvent the doubling problem are the Wilson-type 
and the Kogut-Susskind (or staggered) type discretizations. In the former solution the mass 
of the 15 doublers is increased as compared to the original fermion. This is achieved by adding 
to the naive action Sp^"^ a term that contains a second derivative: —r /2^^%p{n)d^d^-il){n) 
with r an arbitrary constant. This extra term is proportional to a, therefore it vanishes in the 
continuum limit. On the other hand it raises the masses of the unwanted doublers proportional 
to 1/a. The action with Wilson-fermions then has the form 



^Wilson ^Y^^^^j^ Ar)^{n)'4){n) 

n 

{^{n){r - '^^)U^{n)i){n + afi) + ^(n + a/i)(r + lt,)Ul{n)i){n)) 



[2.21) 



This action breaks chiral symmetry for r ^ even for zero quark masses on a lattice with 
finite lattice spacing. This implies that the quark mass will have an additive renormalization 
which makes it very difficult to study chiral symmetry breaking as for that a fine tuning of the 
parameter m is required. 
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Another popular method to get rid of the doublers is to modify the naive action such that 
the Brilluoin zone reduces in effect. This is achieved by distributing the fermionic degrees of 
freedom ipa over the lattice such that the effective lattice spacing for each component is twice 
the original lattice spacing. We lay out the spinor components of il){n) on the sites of the 
hypercube touching the site n. This way, in four dimensions the degrees of freedom reduces 
by 75%, so this formulation describes only 4 flavors of quarks. This discretization is referred 
to as the Kogut-Susskind or staggered fermionic action. 

Applying an appropriate local transformation on the fields if) and ip the naive action (2.19) 
can be diagonalized in the spin-indices a. This way the Dirac-matrices (7^)0/3 are eliminated 
and with the new fields x ^^^d x the staggered fermionic action is 



where the only remnants of the original Dirac structure are the phases Ti^{n) = i)"i+"2-i 
A huge advantage of staggered fermions is that for zero quark mass an U{l)iiX U{1)l symme- 
try (which is a remnant of the full chiral symmetry group) is preserved. Due to this there is 
no additive renormalization in the quark mass and thus no fine tuning - as opposed to Wilson 
fermions - is necessary. Consequently, using the staggered action it is possible to study the 
spontaneous breaking of this remnant symmetry and the corresponding Goldstone-boson. We 
remark furthermore, that the staggered action introduces discretization errors of O(a^). 

As mentioned above, the staggered discretization describes 4 flavors. Since (in the Wilson 
formulation) including a second quark flavor could be realized by inserting another determinant 
in (2.16), one expects that taking the root of detM can be used to decrease the number of 
flavors. Thus for staggered fermions, the following partition function 



is expected to describe Nq flavors. This "rooting" trick [22] is theoretically not well established, 
since the above Z cannot be proven to correspond to a local theory (unlike that in (2.16)). 
Nevertheless numerical results seem to support the validity of the rooting procedure. 

A further possibility to discretize fermions and simultaneously get around the Nielsen- 
Ninomiya theorem resides in defining a lattice version of chiral symmetry: {^,75} = 0(a). This 
way the doublers are avoided and chiral symmetry in the continuum limit is recovered. These 
regularizations are called chiral fermions. Among solutions satisfying lattice chiral symmetry 
are the overlap and the fixed-point fermionic actions; the domain wall fermions on the other 
hand provide an approximation to such a Dirac operator. 



^^tag = - ^ r/^(n) [xin)U^{n)xin + afi) - xin + a/i)f/^(n)x(n)] + ma x(^)x(«) (2-22) 





(2.23) 
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2.3.4 Positivity of the fermion determinant 

In this section an important property of the fermion matrix M = Ip + ml is pointed out. It 
is straightforward to check that any of the above presented fermionic lattice discretizations 
satisfies the condition of 75-hermiticity^: 

= 75^75 (2.24) 

For example for naive fermions (2.19) the Dirac-operator takes the form (with the color and 
Dirac indices suppressed): 



1 4 

IPnm = 2 S {U^{n)-i^5m,n+a(i " Ul{n - ajl)-i f,5m,n~a(i) (2.25) 

^l=l 

Taking into account that for any yU the gamma-matrices satisfy 757^ = —7/^75 and 7I = 1, we 
have 

4 



(75^75)„„ = 2 (-t^M("')7M'5m,n+aA + Ul^U - ttfl)-^ ^6m,n-afi) 



2 

which is indeed the nm matrix element of the adjoint of (2.25), Ipn^- Note that (2.26) is the 
adjoint in color, Dirac, and coordinate space, since the gamma-matrices are self adjoint and n 
and m is interchanged. 

Let us now consider the eigenvalue equation of Ip: 

JPXn = KXn (2.27) 

and define the characteristic polynomial of Ip as -P(A) = det(^ — Al). Then one obtains 

P(A) = det {^l{Ip - Al)) = det {^^{Ip - 



det (jp^ - Al) = det [IP - A*l)* = P*(A*) 



(2.28) 



That is to say, if A is an eigenvalue (P(A)=0), then A* is also an eigenvalue, since P(A*) = also 
holds. This implies that the eigenvalues of Ip are either real, or consist of complex conjugated 
pairs, i.e. the determinant of Ip is real. 



^It is evident that equation (2.24) also holds with Ip replaced by M. 
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As will be discussed in section 2.5, the determinant det M is to be used as a probability 
weight and thus needs to be nonnegative. Combining (2.24) with chiral symmetry (2.20) one 
observes that = —-0, i.e. Ip is antihermitian: its eigenvalues are purely imaginary, A„ = irjn- 
Thus the eigenvalues of the fermion matrix are of the form m ± irjn and the determinant of M 
is indeed a nonnegative real number. 

It is straightforward to prove that (2.24) holds also for (2.21) and (2.22) and as a conse- 
quence the staggered fermion determinant is always nonnegative. Note however, that Wilson- 
fermions do not exhibit chiral symmetry, which means that there can be eigenvalues of M that 
are real and negative which can spoil the positiveness of the Wilson- fermion determinant. 

Also note that as a result of the inclusion of a 6'-term or a chemical potential, the Dirac- 
operator will no more be 75-hermitian. This in turn has serious consequences on the positivity 
of det M, see section 3.3.1. 

2.4 Continuum limit and improved actions 

After the field variables of (2.3) have been discretized the continuum action is obtained by 
carrying out the a — )■ limit. While the discretization procedure is not unique, different 
lattice actions have to give the same continuum limit. Accordingly, the expectation value of 
an arbitrary operator on the lattice can be written as 



with (0) being the expectation value of the operator in the continuum theory and the second 
term is the deviation or 'lattice artefact' caused by the discretization. How fast the discretized 
action converges to the continuum one is determined by the exponent p > (and the coefficient 
of this term). For the Wilson gauge action this scaling is proportional to (i.e. p = 2). For 
an improved action with larger p the scaling is faster. Therefore with an improved action one 
may be able to approach the continuum limit faster, on the other hand, a complicated action 
can significantly slow down the simulation. The optimal choice may depend on the observable 
in question. 

It is easy to see that the continuum limit of the lattice theory is equivalent to a second- 
order critical point of the underlying statistical physical system. Indeed, let us consider a 
particle with a finite mass m. This mass is a physical (constant) number, irrespective of 
how one measures it on the lattice. On the other hand, the mass measured in lattice units 
771 = 7na clearly has to vanish as a — )■ 0, and therefore the corresponding correlation length 




(2.29) 
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has to diverge: ^ — t- oo. This is just the characteristic property of a critical point in statistical 
physics. 

Near the critical point the statistical system exhibits the property of universality. This 
means that in this region the long-range behavior of the system depends only on the number 
of degrees of freedom, the space-time dimension and the symmetries of the theory. Conse- 
quently, the actual form of the action is less and less important; only the relevant operators 
matter. Nevertheless, irrelevant operators (which converge to zero as a — 0) may modify the 
scaling (2.29). 

It was already mentioned that in general, as part of the renormalization program the bare 
parameters of the theory will depend on the regularization. On the lattice this means that 
these parameters will become a function of the lattice spacing, and as one approaches the 
continuum limit, they have to be tuned as a function of a. 



2.4.1 The line of constant physics 



At a fixed temporal size Nt one can change the lattice spacing by varying the bare parameters 
of the action: the inverse gauge coupling /3 and the quark masses rriq. The fact that towards 
the continuum limit the lattice should reproduce the continuum physics, dictates the functional 
relation between these parameters. This relation ensures that for each lattice spacing a "physics 
is the same". A possible way to define this line of constant physics (LCP) is to fix ratios of 
experimentally measurable quantities to their physical value. 

In QCD with 2 + 1 flavors we have three 
independent parameters: /3, rriud and m^. For 
the study of the phase diagram (chapter 4) 
and the QCD equation of state (chapter 5) 
we fix the functions ms{(3) and rriudiP) such 
that the ratios Jk/^tt and Jk/i^k are at 
their experimental value^. Through this pro- 
cedure we get for the ratio of quark masses 
^s/'^ud = 28.15. Note that different defini- 
tions may result in different functions ms{P), 
TTT'udil^), but these differences converge to zero 
as the continuum limit is approached. The 
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Figure 2.1: The line of constant physics. 



^Here fx, ni-K and uik are the kaon decay constant, the pion mass and the kaon mass, respectively, which 
we take from [23]. 
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detailed determination of the line of constant physics can be found in [24, 25]. This definition 
of the LCP was used in the study of the phase diagram; the corresponding ms{f3) and mudiP) 
functions are shown in figure 2.1. For the study of the QCD EoS this relation was further 
improved [12]. 

2.4.2 Scale setting on the lattice 

On the lattice one can only measure dimensionless quantities. A measurement of e.g. the ratio 
of two particle masses can be compared to the physical value of this particular ratio, as it was 
used to define the LCP in the previous subsection. To determine the lattice spacing itself, one 
has to measure an experimentally accessible observable A in lattice units, i.e. in units of a 
certain power d (the mass dimension of the observable in question) of the lattice spacing a: 

^latt = Acxp ■ a'' (2.30) 

Then the lattice spacing can be calculated using the experimental value of the quantity Aexp. 

Arbitrary dimensionful observable can be used to define the lattice scale in this manner. A 
possible choice is one using the static quark- antiquark potential V{r), which can be measured 
on the lattice using spatial-temporal loops constructed from the gauge field (the Wilson loops) . 
In the confined phase the potential is linearly increasing with r. The coefficient of this term is 
given by the string tension a: 

V{r) or (2.31) 

The potential also contains a Coulomb-like repulsion which dominates at small distances. The 
shape of the potential as a function of r can be used to implicitly define an intermediate 
distance r^: 

= = 1.65 (2.32) 

dr / 

The string tension and the parameter ro are only well defined in pure gauge theory. In the 
presence of dynamical quarks, at increasing distances mesons can be created from the vacuum 
and the string between the two color charges can break, making a and tq ill-defined. 

Therefore in dynamical simulations it is more practical to determine the lattice spacing in 
terms of a mass or a decay constant. For the study of the phase diagram (chapter 4) and the 
QCD equation of state (chapter 5) we fixed the scale by measuring the kaon decay constant 

fK. 

A further possibility is to use the critical temperature Tc to fix the scale. To this end one 
has to determine the critical couplings /3c on lattices with various temporal extent Aj. This 
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scheme is particularly advantageous in pure gauge theory, where the phase transition is of first 
order [26-31], and therefore Tc is sharply defined (as opposed to the case of full QCD, where 
a broad crossover separates the phases [4]). Thus in the study of the pure gauge equation of 
state (chapter 6) this approach was followed. 



2.4.3 Symanzik improvement in the gauge sector 



The scaling (2.29) can be improved by inserting further gauge-invariant terms in the lattice 
action. It can be proven that the plaquette is the only relevant operator that can be built 
from purely gauge links. The second simplest combination is the 2x1 rectangle W^^^, i.e. the 
ordered product of links along such a rectangle. The resulting improved action can be written 
as 



oSymanzik 



-/3 



Co J2 ReTrf/;^ (n) + ci J] ReTrf/^^^^l 



n] 



n,fi<u 



(2.33) 



The lattice artefacts of the Wilson gauge action contain O(a^) and (D{g'^a'^) terms. If the 
coefficients in (2.33) are set to cq = 5/3 and ci = —1/12, the O(a^) term is eliminated, thus 
the above combination will approach the continuum theory as 0{a^) on the tree level. This 
action is therefore called the tree-level improved Symanzik gauge action. 



2.4.4 Taste splitting and stout smearing 

As already mentioned, the staggered fermion discretization (2.22) describes (before applying 
the rooting trick) 4 flavors of quarks. The masses of these four fermion species (which are in 
this case called tastes) are however not the same. The SU{4) flavor symmetry is violated by 
the staggered formulation, and as a result each continuum hadron state has a corresponding 
multiplet of states on the lattice: due to the taste symmetry violation the masses of these 
states are split up [32]. This introduces a discretization error which is important mainly at 
low energies [33-35]. 

As an example, 16 lattice states correspond to each continuum pion state, each of them 
contributing with a 1/16 weight. The following table lists the members of the lattice pion 
multiplet with the taste structure (a 4 x 4 complex matrix, Fq,) and the multiplicity {ria)'- 
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Only a = behaves like a Goldstone-boson, i.e. its mass vanishes in the chiral limit. 
The other 15 states have masses of the order of several hundred MeVs for sensible values 
of the lattice spacing. Though these mass differences vanish in the continuum limit, it is 
very important to suppress them as much as possible. The effect of the heavier "pions" on 
thermodynamic observables can be significant: they can reduce the QCD pressure and can also 
shift the transition temperature. 

Strategies for the suppression have been studied extensively. An effective way to reduce 
splitting is to eliminate the ultraviolet noise from the gauge links (which appears as a result 
of the introduction of the finite lattice spacing), and "smear" the links. During the smearing 
process each link is replaced by an appropriately defined average of the surrounding links. One 
possible way is to add to the gauge link the "staples" around it: 



Ufj,{n) Pr 



[n + ah')Ul{n + afi) 



(2.34) 



with p a constant parameter. As a sum of S'f/(3) matrices the result in general will not be 
an element of the gauge group and thus a projection back to SU (3) (denoted above by Pr) is 
necessary. This specific smearing method is called stout smearing [36]. 

The whole process can be repeated several 
times in order to increase the smoothness of 
the links. In the simulations that I present in 
this thesis p = 0.15 was set and the smear- 
ing was carried out twice in a row. Stout 
smearing is proven to significantly reduce the 
lattice artefacts originating from taste split- 
ting. The mass splitting in the pion multiplet 
for the stout smeared action is shown in fig- 
ure 2.2 as a function of the lattice spacing 
(blue lines). For physical quark masses the 
pion state with the lowest mass is adjusted 
to the mass of the continuum pion. For com- 
parison the splitting is also plotted for the 
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Figure 2.2: Masses of lattice pion tastes as func- 
tions of the lattice spacing for the stout (blue lines) 
asqtad improved action [33] (red lines in the and the asqtad action (red lines). The pion states 
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2.4.5 Improved staggered actions 

The O(a^) scaling of the staggered fermionic action can also be improved by considering a more 
complicated discretization for the derivative term in (2.22). Beside the 1-link term^ 

iy^i'0)(n) = ij{n)U^{n)'^{n + afi) (2.35) 

one can also include higher order terms, like all the possible 3-link contributions. These are 
schematically written as the linear (3,0) and the bent (1,2) terms 

Wl^'^\n) = tlj{n)Uf,{n)U^{n + afi)U^{n + 2afi)^{n + 3a/i) 

Wl^^i^\n) = %i){n)U^{n)U^{n + afi)U^{n + afi + ai>)ip{n + afi + 2av) (2.36) 
+ ■ip{n)Uy{n)Uu{n + ai')U^{n + 2au)ip{n + afi + 2au) 

Furthermore, the 1-link terms can also be smeared with a method similar to the one presented 
in subsection 2.4.4. By an appropriate setting of the coefficients of these improvement terms 
(such that the rotational symmetry of the quark propagator is improved [37]) one can achieve a 
better scaling at the tree level (i.e. at zero gauge coupling). This implies that these actions (like 
the p4 or the asqtad action) approach the continuum action faster at very high temperatures. 
On the other hand this improvement does not suppress the taste splitting and therefore large 
lattice artefacts may be expected in the low temperature region (where the lattice spacing is 
large). 

Moreover, the splitting in the tastes can also produce O(a^) errors through taste-exchange 
processes. By a further improvement these processes can also be suppressed; the resulting 
action is called the hisq discretization [38]. The hisq action together with the stout smeared 
action are proven to have significantly smaller splitting between the various tastes as compared 
to the asqtad or the p4 action. At low temperatures these actions are therefore expected to 
produce more reliable results. 

2.5 Monte-Carlo algorithms 

In order to determine the expectation value of some observable (which is necessary to measure 
e.g. the thermodynamic quantities of chapter 3) one needs to calculate the functional integral 
of (2.16). This integral, as discretized on a four dimensional lattice can have a dimension as 
high as 10^, which excludes usual numerical integration techniques. Such integrals can only 

^The staggered fermionic field is denoted here and also in the following by ip. 
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be calculated by Monte-Carlo (MC) methods based on importance sampling. Furthermore, 
because of the Grassmann nature of the fermion fields standard MC methods are not applicable, 
and instead one integrates out the quark degrees of freedom to obtain the determinant of M 
in the partition function, as in (2.16). 

Therefore, the expectation value of an arbitrary observable cf) can be written as 

{^) = ^J DUcp det M e-^°(^) (2.37) 

Importance sampling means that instead of selecting configurations in U space randomly, we 
generate them according to the distribution 

p(U) = ^ det Me-^«(^) (2.38) 

such that we have a set of configurations {t/'-*^} with i = 1 . . . N. Then the expectation value 
of the observable is readily obtained as (assuming that the configurations are independent) 

1 ^ 

(<f^) = :^^ (2.39) 

i=l 

In practice the -> oo limit cannot be carried out since one only has a finite sequence of 
configurations. The deviation from the exact expectation value in this case is given by terms 
of 0{N-^/^) and can be estimated using the jackknife method [19]. 

Note here that in order to interpret p{U) as a probability measure and apply importance 
sampling, the fermion determinant has to be nonnegative. This constraint is fulfilled for the 
staggered lattice action (2.22), if a chemical potential is not present, see section 2.3.4. 



2.5.1 Metropolis-method 

In order to generate configurations according to the desired distribution the only possible way 
is to construct a Markov chain, i.e. to generate the new configuration {U'} from a previous 
one {U} with a probability P{U' ^ U). Markov chains in general converge to the distribution 
p{U) if the above probability fulfills ergodicity (i.e. by successive steps the whole U space can 
be covered) and detailed balance, which means 

p{U)P{U' ^U)= p{U')P{U ^ U') (2.40) 

A simple Markov process is produced by the so-called Metropolis algorithm. Here, first one 
generates a new configuration {U'} by a random change, and then accepts this according to 
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the probability 



PMet(f/' ^ U) 



min 1 , e 



{Sg{U')-Sg{U)) 



det M([/0 
det M(f/) 



(2.41) 



If the new configuration is not accepted, the original configuration remains for the next step^. 
This procedure is however very inefficient, since it involves the calculation of the fermion 
determinant (i.e. {N^NtY fioating-point operations) in each step. Furthermore, the consecutive 
configurations are certainly not independent. 

In this context it is useful to introduce the notion of autocorrelation time, which is the 
number of steps after which the new configuration can be considered independent of the original 
one (i.e. when the correlation of the two falls below some small number). A further important 
quantity is the thermalization time, which can be identified with the number of steps necessary 
for the ensemble of the generated configurations to reach the equilibrium distribution p. 



The Metropolis-algorithm can be improved in many aspects. A much more effective way to 
generate configurations is by means of the so-called Hybrid Monte Carlo (HMC) method [39, 
40], which is a mixture of the Metropolis and the molecular dynamics method. First, we make 
the observation that the determinant of a hermitian matrix H can be written as the (bosonic) 
integral of an exponential: 



The Lp fields are referred to as pseudofermions. The fermion matrix M itself is not hermitian, 
but the combination M'^M obviously is and thus can be used in the above formula in place of 
H to obtain^ 



^Here it can be explicitly seen that the positiveness of the determinant is necessary to obtain a probability 
for which Pe [0,1]. 

^One can further notice that the staggered action (2.22) connects only nearest neighbors and therefore in 
the matrix M only the odd-odd and the even-even elements are nonzero. Furthermore, the determinant 
can be factorized as 

det(A/tA/) = det(MtAf)evcn • det(AfU/)odd 

and using the actual form of the Dirac matrix it is also easy to see that the even and odd factors are equal. 
Therefore (2.43) indeed holds. 



2.5.2 The Hybrid Monte-Carlo method 




(2.42) 




(2.43) 
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Using this representation of the determinant the partition of (2.37) can therefore be written 
in the following form: 

Z=^J D[/Dv?t2)^e-SG(c/)-¥'t(A/tA/)-vU ^2.44) 

In the molecular dynamics method one introduces a new simulation time parameter r and 
considers the time development of the system in this new variable, which can be obtained via 
the Hamiltonian formulation. The canonical variables^ of this Hamiltonian are the gauge fields 
f/„(r) and the corresponding conjugate momenta n„(r), for a fixed value of the pseudofermion 
fields if. Therefore we introduce the conjugate momenta and integrate over them also to rewrite 
the partition function as 



Z = ^ [ 2)n2)t/D(/.t2)^e-ETrn^-5G(J/)-^nAftM), 
CC J 

C'= ! Me-S™" 



(2.45) 



Thus the Hamiltonian of this system can be written as 

^ = ^ E ™„,(r)2 + Scmr)) + ^\M\Un{r))M{Un{r)))~l^^ (2.46) 

n 

The canonical equations of motion as derived from 'K can now be solved as a function of r. 
Along the solutions n„(r), Unir) the "energy" !K of the system is constant. Thus, advancing 
along such a trajectory corresponds to a special Metropolis step for which the acceptance 
probability is 1. In this formulation the expectation value of an observable is obtained by 
averaging along the classical trajectory. 

Since in practice the canonical equations can only be integrated approximately, in some 
discrete steps of 6t, the conservation of energy will also be approximate. A possible prescription 
to carry out this numerical integration is the so-called leapfrog algorithm, which introduces 
errors of 6^ ~ 5r^. However, if a Metropolis acceptance test is inserted at the end of each 
trajectory, the systematic error caused by the finite 6'K can also be eliminated. 

From the numerical point of view, the most demanding part of this algorithm is that in 
each step of the molecular dynamics trajectory (and also in the final Metropolis step), a matrix 
inversion has to be carried out to obtain the momenta H„ (which contains terms of the form 
(M'fM)"V)- Equivalently, one has to exactly solve the system of linear equations 

^ = {M^M)x (2.47) 



Here the index n runs over all the links. 
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This can be solved by e.g. the conjugate gradient method. The time this algorithm needs for 
solving the above equation is proportional to the condition number of the matrix which, in 
turn, is related to the inverse quark mass. Due to this, simulations of systems with smaller 
quark masses are increasingly difficult. 

2.5.3 HMC with staggered fermions 

In order for the staggered lattice action to describe 1 (or 2) quark flavors, one needs to take 
the fourth (or square) root of the fermion determinant, as in (2.23). Unfortunately, in this case 
the conjugate gradient method fails. However, one can use rational functions to approximate 
the root function as 



For each term, the system of equations in (2.47) can now be solved and using J = 10 — 15 terms 
and appropriately tuned coefficients the exact solution for the inversion of the root of M'^M 
can be recovered. This algorithm is referred to as the Rational HMC (RHMC) method [41]. 
This algorithm was used to obtain all of the results presented in this thesis. 




(2.48) 
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Chapter 3 

QCD thermodynamics on the lattice 



After this brief introduction to the lattice approach of QCD, I will analyze the theory from the 
thermodynamic aspect. In this section I will identify the symmetries of the theory and consider 
the corresponding observables that one can use to extract the thermodynamic properties of the 
system. First of all, let us consider the partition function (2.23) in the staggered formulation, 
generalized to the case of a higher number of flavors. The flavors are labeled by g, each having 
a mass of rrig, an assigned chemical potential fig and a degeneracy Ng (thus, in the 2 + 1-flavor 
system A^^ = 2 and = 1). Also, let us denote the total number of quarks as Nq = Ylq-^q- 
After integrating out the quark fields ipg, the partition function reads 



where the dependence of the determinant on the chemical potential and the mass is explicitly 
written out. For each q, the chemical potential fi is treated in the grand canonical approach i.e. 
the action is complemented by a term fiN where A^ is the number of quarks in the system. The 
lattice implementation of the chemical potential is studied in detail in section 3.3. Nevertheless, 
note already here that the chemical potential enters only the fermionic part of the action and 
is not present in Sq- The expectation value of an arbitrary observable (p based on the above 
partition function is written as 



Here and in the following the fermionic determinant is calculated using the staggered dis- 
cretization, which is used for obtaining all of the results in this thesis. 




(3.1) 




DU e-^«(^) Yl M^'^^U, fig, rrig) 



(3.2) 
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3.1 Thermodynamic observables 

In lattice simulations the partition function (3.1) itself is not directly accessible. There are 
on the other hand various observables one can measure using the partial derivatives of log Z. 
Such observables will in general be sensitive to the transition between hadronic matter and the 
QGP and thus play a very important role in thermodynamic studies. These are often referred 
to as approximate order parameters. The reason for this will be discussed in more detail in 
section 4.1. 



3.1.1 Chiral quantities 

It is well known that the QCD Lagrangian (2.3) exhibits an U{Nq)l x U{NQ)ji chiral symmetry 
in the limit rrig — )■ where all flavors are massless. In particular, an axial U{1) transformation 
on any field ipq leaves the Lagrangian invariant: 

V'g ^ e^^^^V'g, i^g ^ tpge''^'^' ^ L{mg = 0) L{mg = 0) (3.3) 

The staggered fermion formulation - although not fully chirally symmetric - is also invariant 
under such a transformation. This part of the group is therefore often referred to as the 
staggered remnant of the full chiral group. Note that Wilson-fermions do not preserve this 
symmetry and thus in this case the breakdown of chiral symmetry would be much more difficult 
to study. 

The order parameter of this symmetry is the quark chiral condensate '(/'gV'ij- The chi- 
rally broken (low temperature) phase is characterized by a nonzero vacuum expectation value 
{'4'q'^q) > 0, while in the symmetric (high temperature) phase {ipqipq) = 0. The chiral con- 
densate for the flavor q can be written as the partial derivative of the partition function with 
respect to the quark mass niq-. 



f'DUe~'oiu)^^,M^./.dlogdetM 



(3.4) 



A Z J dmq 
Using the equality logdet = Trlog and taking into account that = 1, one obtains^ 



i^q^q) = ^ {Tt{M-')) (3.5) 



^Here the identity (M^^)' = -M-^M'M"^ is also used with the prime denoting differentiation with respect 
to an arbitrary variable. Also note that although suppressed, the Dirac operator of course flavor-dependent. 
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The second derivative of log Z with respect to the quark mass is also of interest; it is called 
the chiral susceptibility: 

\Xqq/ 



= _^ (Tr(M-i)>' + ^ (Tr(M-i)Tr(M-i)> - ^ {TiiM-'M-')) ^^"^^ 

where the contribution originating from the single expectation value is divided into a discon- 
nected and a connected part with 



^1rr„n\J-l\ ,,disc. _ fj. „/, ^2 ,,conn. _ ^{ipq'lpq 

Ig 



V^.V'. = ^Tr(M-i), xt^- = (^A)', xr"- = ^^ (3-7) 



Usually one studies the chiral condensate density and the chiral susceptibility density, which 
is obtained from the above combinations after a multiplication with I/V^d = T/V. 

3.1.2 Quark number-related quantities 

A part of the full chiral group is the U{1) vector symmetry. This symmetry of (2.3) corresponds 
to the freedom of redefining the phases of the quark fields 

ijg ^ e^Vg, i^g ^ ^<?e"*' ^ L^L (3.8) 

This U{1) symmetry - which is valid for arbitrary niq - is related to quark number conservation. 
Thus in this regard it is useful to study the quark number density Ug, which is proportional to 
the first derivative of log 2, with respect to the chemical potential /i^: 

_ glogZ 

{n,) = (3.9) 
In the same manner as in (3.4) and (3.5) one obtains: 

K) = ^(Tr(M-iM')> (3.10) 

where the prime indicates a derivative with respect to the chemical potential assigned to the 
quark labeled by q. Second derivatives with respect to the various chemical potentials can 
also be defined. For the study of the phase diagram we will only be interested in the diagonal 
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susceptibilities, i.e. those that are twice differentiated with respect to the same /ig. These 
observables we will refer to as the quark number susceptibilities: 

^2 O- / AT \ , /Ar\2 



(x,>-^-(^) (T.(M-M')>% (^) (MM-Hir) 



^ (Tr(M-^M" - M-^M'M-^M')) 



(3.1i; 



which is compactly written as 

(x,) = -K)' + (xf^- + xr"-> (3-12) 

where a disconnected and a connected term was once again defined: 

N 8n 
n, = ^TT{M-'M'), xf^-=n^, = ^ (3-13) 

To obtain the corresponding densities one should once again multiply by T/V. Moreover, it is 
customary to study the combination (xq) /T^ for reasons discussed in section 4.2.3. 

Note that despite its name, the quark number susceptibility does not exhibit the peak-like 
structure that is usual for a susceptibility in statistical physics. This is due to the fact that Xq 
can also be written as the first derivative (with respect to fig) of the thermodynamic potential 
logX. 

3.1.3 Confinement-related quantities 

In the limit rriq — )■ oo of infinitely heavy quarks the QCD Lagrangian possesses an additional 
symmetry. In this limit quarks decouple from the theory and one is left with a purely gluonic 
system described by Sg- This system is invariant under a center transformation of the temporal 
links, i.e. a transformation where each temporal link Ui{n) is multiplied by a Z G ^(3) 
center element (and the adjoint links by Z'^). Since by definition Z commutes with every link 
variable, any closed loop of links is invariant under this transformation, except for loops that 
wind around the temporal direction. Therefore there is an observable that explicitly breaks 
this Z{?)) symmetry, which can be constructed by multiplying the links along timelines: 

. Nt~l 

P=y Trnf^4(n) (3.14) 

ni,n2,n3 ?14=0 

This observable is called the Polyakov loop. Its expectation value is connected to the free 
energy of a static quark- ant iquark pair taken infinitely far apart Fqq{r — )■ oo): 



34 



3.1 Thermodynamic observables 



In the low temperature phase of pure gauge theory Fqq{r) diverges as r oo and as a conse- 
quence (P) = 0. This is just the phenomenon of confinement: it takes infinitely large energy 
to separate a quark from an antiquark (this can be thought of as the presence of a string that 
connects the quark and the antiquark). At high temperatures quarks are no longer confined 
and thus Fqq{r — oo) is finite, producing a nonzero expectation value for the Polyakov loop. 
In view of the observation that (3.14) is not invariant under ^(3) transformations, this means 
that the Polyakov loop acts as the order parameter of center symmetry: at low temperature 
the symmetry is intact, while at high temperatures it is spontaneously broken. 

In full QCD with finite quark masses 2'(3) is not a valid symmetry anymore, as the fermion 
determinant contains terms that transform nontrivially. This non-invariance can pictorially 
be described by the fact that quark- antiquark pairs can be created from the vacuum and the 
"string" formed between strong charges can break. Nevertheless, the Polyakov loop still signals 
the transition from hadronic matter to the QGP by increasing from almost zero to a larger 
value. 

For the scale setting procedure of chapter 6 we will also use the susceptibility of the Polyakov 
loop, which is defined as 

XP = V{{P')-{Pf) (3.16) 
3.1.4 Equation of state- related quantities 

The partition function also serves to define observables that can be used to establish the 
equation of state of the theory. Such observables play an important role in describing the ther- 
modynamic properties of the system; their definition is given in this section. These definitions 
will be applied in chapter 5 for the determination of the equation of state both in pure gauge 
theory and in full QCD. 

The free energy density is related to the logarithm of the partition function as 

/ = -^logX (3.17) 

The pressure is given by the derivative of T log 2. with respect to the volume. Assuming that we 
have a large, homogeneous system, differentiation with respect to V is equivalent to dividing 
by the volume. Therefore in the thermodynamic limit the pressure can be written as minus 
the free energy density: 

p = - lim /. (3.18) 

V^oo 

In lattice simulations the validity of this assumption has to be checked. This will also be 
elaborated on in chapter 5 and chapter 6. Having calculated the pressure as a function of the 
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temperature p{T), all other thermodynamic observables can also be reconstructed. The trace 
anomaly J is a straightforward derivative of the normalized pressure: 

I^e-3p = T'^P-^ (3.19) 

This combination is often called interaction measure as it measures the deviation from the 
equation of state of an ideal gas e = 3p. The inverse relation can easily be written as 



piT) } I{T') 



2^4 / 





dT' (3.20) 



Using the pressure and the trace anomaly the energy density e, the entropy density s and the 
speed of sound Cg can be calculated as 



7 + 3^ s-i±^ c^-^ 



(3.21) 



Note that in the absence of a chemical potential all the above observable are functions of only 
one variable, namely the temperature T. Therefore varying the pressure or the energy density 
inevitably modifies the entropy density also. Nevertheless, the speed of sound remains a well- 
defined quantity, since it can be rewritten as a ratio of two partial derivatives at constant 
volume [42]: 

2 dp 



' dT 

3.2 Renormalization 



df 



(3.22) 



In the previous sections I presented the most important observables that one can use to study 
the thermodynamic properties of the system. These are all constructed using the free energy 
density, which is however an ultraviolet divergent quantity. Thus, in order to have a meaningful 
continuum limit, a proper renormalization to these observables has to be applied. Based on 
dimensional reasoning, the free energy itself contains additive divergences in the following 
form [43]: 

f = fr + 0{a-^) + 0{m\-^) + 0(m^ log(a)) (3.23) 

with fr being the renormalized free energy density. Note that a term T'^a~'^ - despite having 
the correct dimension - is not present in the above expression. The absence of this contribution 
expresses the fact that divergences are in general independent of the temperature. Once one 
calculated the divergences at T = and carried out a proper renormalization, the effect of 
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heating up the system to a finite temperature is just equivalent to assigning different weights 
to the states according to the Boltzmann-factors. This of course should not bring in new 
divergences. 

One can eliminate the additive divergences by subtracting the T = contribution from the 
free energy. Thus the following expression is ultraviolet finite: 

T) = T) - f{ai, T = 0) (3.24) 

where ttj are the bare parameters of the theory, on which the dependence of / is explicitly 
written out in order to emphasize that the subtraction has to be carried out at the same 
value for each bare parameter. Remember that the temperature can be set according to the 
number of lattice sites A^^ in the temporal extent, see (2.11). Therefore, on a finite lattice, zero 
temperature can never be realized. Usually, however, a large enough Nf (such as Nt > A^^) 
represents a system that is so deep in the low-temperature phase that can it can safely be 
considered as having effectively zero T. Accordingly, in the following, T = should be thought 
of in this sense. 



3.2.1 Renormalization of the approximate order parameters 

The chiral condensate - as a result of being a derivative with respect to the bare mass - also 
contains a multiplicative divergence besides the additive one. This can be cancelled if one 
multiplies with the bare mass such that the resulting combination is a renormalization group 
invariant. In order to obtain a dimensionless expression a simple normalization by some mass 
scale can be carried out^: 



— Tn — — 

^MT) = ^ ■ {^^{T) - ^^{T = 0)) (3.25) 



in the same manner, the chiral susceptibility is renormalized with the square of m: 



2 

m 



X^^,r(^) = gi ■ (x^^T) - x^^{T = 0)) (3.26) 

For the normalization one can use the fourth power of the T = pion mass m^, or any other 
(possibly temperature dependent) combination, e.g. or m^T^. This will be elaborated on 
in more detail in section 4.2.3. We remark here that this renormalization procedure leads to a 



^The chiral condensate ipg^'g is abbreviated here as tpijj; the renormahzation procedure apphes to any of 
the quark flavors. Also, dependence on the bare parameters is suppressed from now on. 
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somewhat unusual chiral condensate which vanishes at T = and reaches a negative value at 
T ^oo. 

The quark number density and the quark number susceptibility inherit no divergent con- 
tributions from the free energy density and thus remain finite as the continuum limit is ap- 
proached. The origin of this is the absence of a yU^a~^ term in (3.23), which follows from the 
fact that quark number is a conserved quantity and thus needs no renormalization. Conserved 
currents associated with non-Abelian symmetries are in general not renormalized since they 
obey nonlinear commutation relations and thus their overall normalization is already fixed. For 
Abelian symmetries like the U{1) symmetry of quark number conservation nonrenormalization 
can be deduced from the corresponding Ward identity. 

The divergences of / also show up in the Polyakov loop. These can be eliminated by the 
renormalization of the static quark- ant iquark potential V{r) at T = [44], namely, prescribing 
the condition K(^o) = for the renormalized potential. This can be interpreted for the 
Polyakov loop to satisfy the renormalization condition 

|(P^)| = |(p)|e^(^^«)/2^ (3.27) 

The potential V{r) on the other hand can be determined at T = from Wilson loops. 

3.2.2 Renormalization of the pressure 

On the lattice the renormalization (3.24) can be written in the form 

fl^\a,, iV„ N,) = f^\a,, iV„ N,) - f^\a,, iV„ iV-^) (3.28) 

where the temperature (to which the left hand side corresponds to) is determined through T = 
{Nta{ai))~^ with the lattice spacing treated as a function of the bare parameters. Furthermore 
it is useful to define the integer parameter ^ > 1 by 

^sub _ ^j^^ ^3 29) 

that is to say, the lattice used for the subtraction corresponds to a temperature of T/^. Instead 
of using a large ^ (i.e. a large Nf^^, which is computationally very demanding), one can 
also calculate the renormalized observables using a reasonably small ^ (due to the fact that 
divergences are in general independent of the temperature). This approach is presented here for 
the case of the pressure, but it can also be generalized to any other quantity that is additively 
renormalized. 
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Let us introduce the following combination [15]: 



Pi{T)=p{T)-p{T/0 



(3.30) 



with ^ being an integer larger than one. Using this difference the renormalized pressure can 
be built up as 

p,(T) = p{T) - p{T/0 + p(T/0 - p{T/e) + PiT/e) - p{T/e) + ■■■ 
= p^(T)+p^iT/0+pdT/e) + ... 
Furthermore, for the dimensionless combination this means 

Pr{T) 



(3.3i; 



rp4 



EL 

2^4 



1 



Pi_ 

2^4 



+ 7 



El 

ji4 



+ 



(3.32) 



Since the forthcoming terms are suppressed by increasing powers of (and since the pressure 
itself gives smaller contributions at smaller temperatures) in practice one can safely truncate 
the series after a few terms. Simply using the first term in the series with ^ = 3 causes a 
relative error of somewhat more than one percent. For the study of the EoS in section 5.3 this 
ratio is below the typical statistical error and thus can safely be ignored at the level of the 
present statistics. 

The renormalization prescription (3.32) can also be applied to the trace anomaly /, the 
energy density e and the entropy density s. However, since these can all be derived from the 
pressure, if once the renormalization of the pressure is carried out, the other observables are 
also going to be finite in the continuum limit. 



3.3 Chemical potential on the lattice 

In heavy ion collisions usually there is a nonzero net quark density in the system resulting 
from the initial excess of quarks over antiquarks. In the grand canonical approach to statis- 
tical physics the density of quarks can be controlled by a chemical potential fi, which can be 
introduced by including in the action a term fiN with being the number of quarks. In the 
Euclidean formulation the equivalent of is the four- volume integral of = ip'y^ip. A naive 
inclusion of the chemical potential by an additional /ij4 term in the action would however cause 
quadratic divergences in the energy density [45]. 

However, the term ipig'-jyAyip in the QCD Lagrangian shows that the chemical potential 
actually acts like the fourth component of a purely imaginary, constant vector potential (/i ~ 
igAis). Consequently, a straightforward way to introduce n on the lattice is to complement the 



39 



3. QCD THERMODYNAMICS ON THE LATTICE 



fourth component of the gauge field with it in the fermionic part of the action^. This amounts 
to multiplying the links in the forward Euclidean time direction with e"'*, and the backward 
links with e~"''^. As a result the fermionic action will be (using the staggered formulation): 



sta 1 r ^ - 



2 

+ r/4(n) [%l}{n)Ui{n)e''^%l){n + al) - %l){n + al)Ul{n)e-''^i){n] 
+ ma ijj{n)^{'i 



(3.33) 



in 



On the other hand it is easy to prove that the fermion determinant can be written as a sum 
over closed loops. In such loops factors of e^"'^ once again cancel, unless the loop winds around 
the Euclidean time direction. For a loop with number of windings the total contribution 
will therefore be (e='='^'^)'^™^'. This however implies that implementing the chemical potential on 
the lattice can be done by multiplying the time-like links by e'^'^'^* = e^^'^ on a single timeslice 
(and the links in the opposite direction in that timeslice by e'^^"^). This implementation 
is connected to that in (3.33) by a transformation of the links, which leaves the fermion 
determinant unchanged. 

The action (3.33) describes one quark flavor. For the case of more flavors a different chemical 
potential has to be assigned to each quark field. As mentioned before, in the energy range 
under study it is useful to take into account the three lightest flavors u, d and s. In heavy ion 
collisions the strangeness of the initial states is zero, and - due to the strangeness-conserving 
nature of the strong interactions - strange quarks s can only be produced together with their 
antiquarks s. This implies that the net strange density throughout the whole process is zero, 
thus /i^ = 0. Furthermore, it is also realistic to set the chemical potential assigned to u and 
d the same: fiu = fJ'd = fJ'Q- This quark chemical potential is then equal to one third of the 
baryonic chemical potential: /ig = /iB/3. In the following, if not stated otherwise, /i will 
always denote the chemical potential assigned to the light quarks /iq. 



3.3.1 The sign problem 

In section 2.3.4 it was shown that the fermion matrix satisfies the condition (2.24). However, if a 
chemical potential is present, the important property of 75-hermiticity is lost, and consequently 

^Remember that the gauge field is defined as = e*'"^^" . 
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the determinant of M will be complex. This implies - as can be explicitly checked using the 
discretization (3.33) - that in such cases the fermion matrix satisfies (for any G C): 



This of course gives the original condition for /i = 0. Note however, that 75-hermiticity also 
holds for purely imaginary values of the chemical potential. 

For a real chemical potential on the other hand the determinant will be some complex 
function and thus cannot be used as a probability weight in (2.16). Since physical observables 
have real expectation values, it is instructive to use the real part Re det(M)e~'^'^ as the weight. 
This can however still be negative and as a consequence can cause large cancellations when 
averaged over different configurations. This is called the sign problem. Note that a nonzero 
imaginary part of the determinant is necessary in order to describe a system with nonzero 
baryon number [46]. We remark furthermore, that the sign problem is a general characteristic 
of Monte-Carlo studies of fermionic systems and is not particular to the lattice approach. 

Recently several methods were developed to circumvent the sign problem and thus access 
the region of small chemical potentials. They are all based on simulations at zero or purely 
imaginary chemical potentials where as was argued for according to (3.34) the sign problem is 
absent. 

In the reweighting method one generates configurations with the action with zero (or purely 
imaginary) /i, and then assigns new weights to them in a way that they describe a /i > 
ensemble [8, 47-50]. Although this approach is exact in the infinite statistics limit, on a 
finite number of configurations the weights oscillate strongly, resulting in a large cancellations. 
Furthermore, since it requires the evaluation of the fermionic determinant, the reweighting 
method is unfortunately restrained to rather small lattices. Another method is to carry out 
measurements at various values of the purely imaginary fi and then analytically continue to 
a real chemical potential according to some ansatz function [51-60]. One needs to carefully 
choose the actual form of the ansatz function since the continued results depend very strongly 
on it. A further approach is the use of the canonical ensemble, where one works in sectors with 
fixed baryon number [61-63] and once again the determinant of the fermion matrix needs to 
be calculated. Finally, a possible method to analyze the system with nonzero /i is to carry out 
a Taylor expansion around zero (or purely imaginary) chemical potential. Such studies can be 
improved by measuring higher order coefficients in the expansion. This approach can be shown 
to be just the expanded version of the reweighting method. Furthermore, the Taylor-expansion 
is not restricted to small lattices, which allows for the systematic study of finite size and lattice 




(3.34) 
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discretization errors. In principle the expansion is expected to converge up to the singularity 
on the complex /i-plane that is closest to the origin. 

In this thesis I present results regarding the phase diagram that were obtained using the 
Taylor expansion technique. Therefore the next section is devoted to the detailed study of this 
approach. 



3.3.2 Taylor expansion in fi 

It is instructive to begin with two general remarks. First, time reversal symmetry ensures that 
Z{fi) = Z{—fi) i.e. the partition function - and thus e.g. the chiral susceptibility also - is an 
even function of the chemical potential. As a consequence of this, derivatives with respect to 
/i^ and second derivatives with respect to n are proportional to each other at zero chemical 
potential: 



{Xi,i,) 



2 ^i^^^) 



fi=0 



(3.35) 

M=0 



Second, let us expand the fermion determinant in the chemical potential as 

det M(/i) = ao + Oi/i + a2/i^ + as/i'^ + . . . (3.36) 

Equation (3.34) implies that for the determinant det M(/i) = det M^(—/i*) holds. For the 
coefficients G C this accounts to 

ao = Qq] ai = —al; a2 = 0*2] ■■■ (3.37) 

which means that even coefficients a2n are real, while odd coefficients a2n+i are purely imaginary 
(for any ri G N). 

By means of the Taylor expansion approach one reconstructs observables at finite chemical 
potential using their derivatives at zero chemical potential. Instead of yU = 0, the expansion 
could also be carried out around a purely imaginary yU = ZyU/. However, for the curvature of 
the pseudocritical line, which will be the primary objective of the next chapter of the thesis, 
the behavior of observables around /i = is of interest. 

Suppose now that we have an observable 0, which may depend explicitly on /i. Remember 
that /i now denotes the light chemical potential. As a straightforward calculation shows, the 
derivative of the expectation value (0) (see (3.2)) with respect to /x can be written as: 

^ = - (n) (0) + (0') (3.38) 
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where n is the light quark number density as in (3.10) and the prime stands for the derivative 
with respect to the light chemical potential. The second derivative is similarly given by: 



dp? 



= (0(^2 + n')) - (n) {n(p) - (0) {n^ + n') - {n(p) (n) + 2 {nf (0) 
+ (n0') - (0') (n) + (0") + {n<P') - (0') {n) 



(3.39) 



where the combination n^ + n' = T^disc._|_^conn. jg ^.j^^ disconnected plus the connected part of the 
light quark number susceptibility (see definition in (3.12)). Now we can take into account that 
according to (3.37), odd derivatives of det M - such as n - are purely imaginary. Therefore, 
expectation values of such operators necessarily vanish at zero /x. Using this fact we can rewrite 
the above two expressions as: 

^^'^^ = {n<P) + (0') (3.40) 

fi=0 

= (0(n2 + n')) - (0) + n') + 2 (r20') + (0") (3.41) 



dHct>) 



fi=0 



For observables that do not depend explicitly on the (light) chemical potential, like (j) = P 
or = rij, (the latter only depends on fig), the last two terms in (3.41) vanish. However, for 
/x-dependent observables, like = ipip = "^i^ud this explicit /i-dependence shows up and we 
have a nonzero contribution. 

For the study of the phase diagram in chapter 4 two observables will be addressed: the 
renormalized chiral condensate and the strange quark number susceptibility. For = ipip^ the 
formula (3.41) is directly applicable. Note that the additive renormalization (see section 3.2) 
of the condensate does not influence the derivative in question. For the strange quark number 
susceptibility (as in (3.12) with q = s) - even though at zero chemical potential (n^) = - 
there is an additional term that contributes to the second /i-derivative. Finally, at /i = one 
obtains for both observables 

T d'^{%l)%l)r) 

V 

V 



1-1=0 



^ [{tPMn^ + n)) - (#r> {n^ + n')+2 (mPtPl) + (3.42) 



T 1 



[{Xsin' + n'))- ixs) {n' + n') - 2 



^l=o 



(3.43) 
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Chapter 4 

The phase diagram 



As it was pointed out in the introduction, the phase diagram of QCD plays a very important 
role in high energy physics, both from the theoretical and the experimental point of view. The 
diagram displays the phases of strongly interacting matter as a function of the state parameters. 
As argued for in section 3.3, the density of quarks is controlled by a chemical potential and 
the relevant parameter (in e.g. heavy ion collisions) is the light baryonic chemical potential 
fiB- Consequently, it is instructive to study the phase diagram on the hb — T plane. In the 
following the light baryonic chemical potential is denoted simply by 

For yU = various lattice results can be found in the literature. In particular the nature 
and the temperature of the transition were studied in detail. Before I present results regarding 
the phase diagram it is necessary to get acquainted with these = studies, therefore an 
overview of the topic is presented in section 4.1. Based on this grounding next I will address 
the case of nonvanishing chemical potentials. Unfortunately, here the situation is much more 
difficult, since the theory suffers from a sign problem, see section 3.3.1. In section 4.2 I report 
results on the curvature of the transition line, as the first /i > study in which a continuum 
extrapolation was carried out. 

For this study we used the Symanzik improved gauge action (see section 2.4.3), while the 
fermions are discretized using the one-link staggered action with stout-smeared gauge links 
(see section 2.4.4). 

4.1 The finite temperature transition 

For the study of the finite temperature transition the first task is to identify the transition 
itself. To this end one has to explore the transition region and measure some observable that 
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is sensitive to the transition. The singularity or analytic jump in such observables indicates 
that the system transforms from one phase to another. We know that there are two distinctive 
symmetries that are connected to the finite temperature transition: chiral symmetry and Z{3) 
symmetry (see detailed discussion in section 3.1.1 and 3.1.3). The breakdown of the Z{3) 
symmetry and the restoration of chiral symmetry are the two phenomena that signal the finite 
temperature transition. On the other hand, in nature these symmetries are explicitly broken. 
While chiral symmetry is only valid when the quark masses are zero, Z{3) symmetry is present 
only in the absence of dynamical quarks, i.e. when the quarks are infinitely heavy. 

One sees that the transition depends very strongly on the actual values of the quark masses. 
In nature quarks are neither massless nor infinitely heavy, thus the two observables related 
to the two symmetries - the chiral condensate and the Polyakov loop - are not real order 
parameters. Nevertheless, they still signal the transition from one phase to the other and thus 
are useful to study as being approximate order parameters. 

It is instructive to analyze the transition from the 
point of the quark masses a bit more in detail. One 
may can summarize the dependence of the nature of 
the transition on the light and strange masses in a 
compact form on the so-called Columbia-plot (see fig- 
ure 4.1). In the upper right corner, which corresponds 
to — > oo, rriud — > oo, lattice results unambiguously 
identified a first-order transition [26-31]. For the oppo- 
site limit one is restricted to model calculations. Such 
models have the same symmetries as QCD, and uni- 
versality guarantees that near the critical points their 
behavior essentially coincides with that of QCD. At 
the origin (m,, = rUud = 0) we expect to have a first- 
order transition, while along the vertical axis a second- 
order transition is predicted, with a universality class 
of 0(4). The physical point rn^'"^ , m^^^'^ on the other 
hand lies somewhere in the crossover region, as shown by lattice calculations. There are also 
indications that the second-order line in the lower left corner of the plot belongs to the Z(2) 
universality class [64]. The position of this transition line relative to the physical point is also 
of phenomenological interest [9, 55]. Lattice results predict that the crossover transition may 
persist down to ~ O.lm^^^^ [64] (pointed to by the arrow on the plot). Note that the scales on 




mud 

Figure 4.1: The Columbia-plot: the na- 
ture of the finite temperature transition 
as a function of m^ci and ms . First-order 
transitions (green, dashed regions) are 
separated from crossovers (blue region) 
by second-order transitions (solid lines). 
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the plot are exaggerated. It is useful to remark here that the strength of the transition is a 
non-monotonic function of the quark masses. A similar non-monotonic behavior as a function 
of the chemical potential is therefore not to be excluded. 

In the following I concentrate on the physical point, that is to say, the quark masses 
are set to their physical values (see section 2.4.1). The observables of interest are shown 
as a function of the temperature in figure 4.2. From left to right the renormalized chiral 
condensate, Polyakov loop, and the normalized strange quark number susceptibility are plotted 
in the transition region, measured on Nt = 10 lattices. While ipipr decreases from zero to a 
negative value (the details of the renormalization of the condensate is discussed in section 3.2), 
P and Xs/T"^ increase from almost zero to a nonzero value. Apparently there is no singularity 
in the observables shown in figure 4.2. This is no surprise, since on a finite lattice no real 
phase transition can be realized, due to the fact that the partition function 2. is the finite 
dimensional integral of a positive integrand. A real phase transition can only be recovered in 
the infinite volume limit V ^ oo. It is therefore instructive to carry out measurements on 
various lattices with increasing (three dimensional) volume, and study the finite size scaling of 
the observables. Such a careful analysis was carried out recently and the results show that the 
continuum extrapolated observables show no singularities even in the infinite volume limit [4]. 
This reflects that the transition is an analytic crossover, as opposed to being a first-order phase 
transition. 
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Figure 4.2: The renormalized chiral condensate, Polyakov loop and the strange quark number 
susceptibility as functions of the temperature at /x = 0. 



Such a crossover transition - just like the melting of butter - does not exhibit a unique 
critical temperature. On the contrary, one can characterize the finite temperature transition 
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of QCD by a pseudocritical temperature and a (nonzero) width^. In particular, the transition 
temperature and the width can be determined by means of e.g. a cubic fit: the former is related 
to the inflection point of the curve, while the latter to the inverse slope at the inflection point. 
Note that the Polyakov loop and the strange susceptibility tend to give a somewhat higher 
transition temperature as compared to the chiral condensate (see figure 4.2). This difference 
also persists in results extrapolated to the continuum limit [35, 65]. This behavior is not 
unexpected and its reason lies in the crossover nature of the transition. Different definitions 
may give different results for the transition temperature, simply due to the nonsingular behavior 
of the logarithm of the partition function. 

Throughout this chapter the stout smeared staggered fermionic action is used. This dis- 
cretization was applied in the above cited studies also. The transition temperature was however 
also studied using various other formulations, including improved staggered actions or further 
fermionic implementations, see e.g. Refs. [66, 67]. 

4.2 The curvature of the Tc(/i) Hne 

While at zero chemical potential lattice calculations provide reliable and accurate results [24, 
25, 68, 69], there are various questions about the fi > region of the phase diagram still to 
be answered. Two possible scenarios are illustrated in figure 4.3. As discussed above, the 
transition at /i = is a crossover [4] and the transition temperature Tc is expected to decrease 
as we increase /i (one might think that a system with high density favors the deconfined 
phase already ay smaller temperatures). Apart from the actual shape of the Tc(/i) function, 
a particularly interesting question is whether the transition becomes weaker or stronger as fi 
grows. 

A strengthening of the transition could lead to the existence of a critical point, where the 
crossover transforms into a true phase transition (see left side of figure 4.3). Indications for 
the existence of such a critical endpoint are present in the literature, see e.g. [8]. 

Another possibility is that the transition weakens with increasing /i (see right side of fig- 
ure 4.3). Recently it was found on relatively coarse lattices that the transition at small fi is 
weaker as compared to the fi = situation [55, 59]. The existence of the critical point is not 
ruled out by this result, but it requires a non-monotonic behavior. Note, that the strength of 
the transition is also non-monotonic as a function of quark masses as it was demonstrated in 

^In the following the term transition temperature will refer to the pseudocritical temperature. 
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Baryonic chemical potential Baryonic chemical potential 



Figure 4.3: Two possible scenarios for the QCD pliase diagram on the /U — T plane. The left panel 
shows a phase diagram with a transition growing stronger and possibly even turning into a real, 
first-order phase transition at a critical endpoint. The right panel on the other hand corresponds to 
a scenario with a weakening transition and no critical endpoint. The paths representing systems that 
describe the early Universe and a heavy ion collision are also shown by the arrows. 



section 4.1. A similar non-monotonic behavior can be observed at finite /i in a linear a model 
with two quarks [70]. 

The transition temperature Tc{fi) is determined in the following as a function of the chemical 
potential through a Taylor-expansion technique. A particularly important property of this 
function is its curvature n - i.e. its second derivative with respect to /i at yU = 0. Due 
to the crossover nature of the transition, different definitions may give different results for 
Tc(yu) - and thus also for n. To extract this quantity I will show the analysis of two different 
observables that are sensitive to the transition: the chiral condensate and the strange quark 
number susceptibility. Furthermore, the method also reveals information about the change 
in the slope of these observables as a function of T near the transition temperature as /i is 
increased. This is closely related to the change in the width - or, in other words, the strength 
- of the transition. This information can indicate whether the transition grows and turns into 
a real, first-order phase transition or remains an analytic crossover. 



4.2.1 The transition temperature at nonzero /i 

Due to the fact that the partition function Z is an even function of the chemical potential, the 
transition temperature also depends only on fi"^. Accordingly, let us parameterize the transition 
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line in the vicinity of the vertical /i = cLXlS clS 

T,(/i2) = (1 - K ■ ^?/Tl) (4.1 
with Tc being short for Tc(0). This implies that the curvature can be written as 



K 



d(/i2 



(4.2) 

/i=0 



In order to determine the curvature we have to measure Tc as a function of /i for small chemical 
potentials. A straightforward way to do this would be to calculate the value of some observable 
as a function of T at a nonzero A/i using its value at /i = and the derivative (3.41). Then 
one can analyze the behavior of this observable in the transition region by means of e.g. a 
cubic fit to find the inflection point (like shown in section 4.1). However, this approach turns 
out to be somewhat ineffective since then a fitting of the reweighted data is required. Instead 
of this, we use a definition of T^. which is more suitable for determining the curvature. 

Let us consider a quantity 0(T, /i^) that is monotonic in T in the transition region, and 
fulfills the following constraints^: 

hm T^0(r, /i^) = 0, hm T^0(T, /i^) = (4.3) 

that is to say, does not depend on the chemical potential in the limiting cases T — )■ and 
T — )> oo. For any fixed /i we can implicitly define a transition temperature Tclfi"^) as the 
temperature at which (j)(T,fi'^) takes the predefined constant value C: 

<^(^'/^')It=W)=C'- (4.4) 
We will choose a C that corresponds to the inflection point of 0(T, 0)^. 

4.2.2 Definition of the curvature 

Now let us determine the curvature using the definition of Tc(/i^). The total derivative of the 
observable 0(T, /x^) may be written as 



dT 



At=0 



9(/i2 



■ d/i^ (4.5) 

At=0 



^In order not to complicate the notation, in the text wc simply write 4> instead of (0). 
^The dependence of on the constant C is suppressed in the following. 
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Along the Tc{fi^) line, is constant by definition, thus d0 = 0. One obtains 



d/i2 



fj.=0 



dT 



/i=0 



(4.6) 
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Thus, for every C we can define a curvature. Since the Tc(C) function is invertible for the 
whole C range, we can also write the right hand side of (4.6) as a function of the temperature, 
R{T). 

The function -R(T) is related to the distance 
that the 0(T) curve shifts along the T axis as the 
chemical potential is varied. Given 0(T) and -R(T) 
at zero chemical potential, the shift for non-zero /i 
at leading order is -R(T) ■ ^ (the curve moves to 
the left if -R(T) is negative and to the right oth- 
erwise). This behavior is illustrated in figure 4.4. 
Using -R(T) we can define a temperature dependent 
curvature according to (4.2) as k(T) = —T^. ■ R{T). 
The meaning of n{T) is again simple: it gives the 
curvature of the <} 
T at /i = 0. 

We use the value of k(T) at T = Tc to define 
the curvature for a given observable. The shape of 
the k(T) function also has important consequences. 
The slope of k{T) around Tc is related to the width of the transition as follows: if the slope 
is zero, i.e. k(T) is constant around Tc, then all points shift the same amount along the T 
axis when a small chemical potential is switched on. This means that to leading order in /i 
the shape of the (f){T) function (and thus, the width of the transition) does not change. If the 
slope is positive, then points with larger T shift more than the ones with smaller T resulting 
in a compression of points, i.e. a narrower transition. Similarly, a negative slope indicates a 
broadening of the transition. 

Putting all this together, the expression dn/dT is related to the relative change in the 
width W{fi) of the transition as the chemical potential increases: 



T (MeV) 

const, curve which starts from Figure 4.4: Illustration of the behavior of 

the observable (j) at fj, = and fj, > 0. The 
function Tc{iJ?') is defined as the tempera- 
ture where (j){T,fi'^) crosses a constant value 
C (see definition in text). 
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where we assume that W is proportional to the inverse slope of the quantity in question: 

W ^\{d<l>/dT\^^^J-'\. 

4.2.3 The /i-dependence of the measured observables 

The two quantities we used to play the role of (p are the renormalized chiral condensate and the 
strange quark number susceptibility. The definition and renormalization of these observables 
is detailed in section 3.1. For the case of the strange susceptibility it is useful to study the 
combination Xs/T^, since it obeys the conditions of (4.3). It is easy to see that at T = one 
gets Xs/T"^ = and at T — oo the normalized quark number susceptibility Xs/T"^ reaches its 
fi independent Stefan-Boltzmann limit of 1. 

The situation with the chiral condensate is more delicate. The final normalization of this 
observable was carried out by = m'^ (see renormalization condition (3.25)). Note that a 
division by which would also render the condensate dimensionless, changes the temperature- 
dependence to be non-monotonic, which would be disadvantageous in the present context. 
Furthermore, as it was already pointed out in section 3.2, ipipr is zero at T = and approaches 
a negative value as the temperature is increased. A more conventional condensate which is 
positive at T = and goes to zero at large temperatures can be obtained by a constant shift 
which is irrelevant for our present study. 

Finally we need to show that the renormalized chiral condensate is independent of /i at 
T = and T — )■ oo. At T = the partition function is independent of /i as long as fi is smaller 
than a fic critical value (the approximate baryon mass) and no baryons can be created from 
the vacuum. Only for fi > fi^ does the partition function have a non-trivial //-dependence. 
Therefore, for fi < fi^ all derivatives of Z (thus ipipr) are independent of fi. The chemical 
potential regime of interest lies in this region. In the Stefan-Boltzmann limit (T — )■ oo) the 
//-independence can be proven by a straightforward calculation. 

Now, according to (4.1) and (4.6), in order to measure the curvature two derivatives are 
necessary. The derivative d {(p) / dT is calculated numerically, using the /i = data as a 
function of the temperature. The combination 9(0) /9(/t^) on the other hand is determined 
using the technique detailed in section 3.3.2, in particular, see equations (3.42) and (3.43). 
The derivatives of the explicit dependence of the chiral condensate on /i (see last two terms 
of (3.42)) were calculated numerically, using a purely imaginary chemical potential A//j. The 
value of A/tj was varied in the range 0.01 . . . 0.0005, and it was checked that the finite differences 
converge fast enough to the A/ij — > values and the error coming from this approximation is 
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negligible compared to statistical errors. Taking into account these considerations A/ij 
was used. 



0.001 



4.2.4 Finite size effects 
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In order to check for finite size effects the mea- 
surement of the derivatives d {(f)) /d{fi'^) was also 
performed on lattices with different physical 
spatial extent for Nf = 6. In particular we com- 
pared results at /3 = 3.555 obtained on 24^ x 6 
and on 18'^ x 6 lattices. This value of /3 cor- 
responds to about 155 MeV, i.e. is near the 
pseudocritical temperature. The larger box is of 
physical size ~ 5 fm. We observe a good agree- 
ment as the results for d(t)/d{ix^) agree within 
statistical errors for both the renormalized chi- 
ral condensate = ipifjr and the strange quark 
number susceptibility = Xs/T"^^ as can be seen 
in figure 4.5. Thus we conclude that finite size N 
errors can be neglected at the present statistical 

accuracy. Note that this volume-independence is a consequence of the crossover nature of the 
transition. 
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Figure 4.5: The derivative of the observables 
with respect to fi'^ measured on A'j = 6 lattices 
24 (blue) and = 18 (red). 



4.3 Results 

Since the actual shape of the k(T) function is unknown it is instructive to carry out a Taylor 
expansion around Tc in the t = (T — Tc)/Tc dimensionless variable: 

fi:(T) = k(TJ + Co-t + Ci-f (4.7) 

For each lattice spacing (i.e. each Nt) we have several simulation points, corresponding to 
different temperatures. In order to fit all of our points at once, we allow a lattice spacing- 
dependence for the constant and linear terms (having a lattice spacing-dependence of the 
quadratic term is also possible, but it does not improve the quality of the fits). Therefore we 
fit all of our simulation points with the following function: 

k{T; Nt) = k{T^; cont) + co ■ t + ■ f + C2/N^ + Cg ■ t/Nf (4.8) 
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with fit parameters k(Tc; cont), cq, ci, C2 and C3. The independent data points as well as the 
fitted curves (for each Nf and in the continuum) are shown in figure 4.6. 
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Figure 4.6: The curvature k{T) determined using the strange quark number susceptibility (left) and 
the renormalized chiral condensate (right), respectively. A result of the combined fit (see detailed 
description in text) is shown by the gray band. The fit results for the individual Nt = 6, 8, 10 lattices 
are shown by the red, blue and green curves, respectively. 



The x^/d-o.f. values of the two fits are 1.08 and 1.29, respectively, indicating good fit 
qualities. The continuum curvatures are given by the k{Tc; cont) fit parameter. The final 
results concerning the curvature and the relative change in the width of the transition are 

^(Xs/T2) ^ g_gggg(^;L4)^ ^{^^r) ^ 0.0065(20) 

AW^/H^fc/^') = 0.033(16), AW/W^^^^^ = 0.030(18) 

with the statistical error of the appropriate quantities in parentheses. The curvatures obtained 
from the two observables are consistent with each other within errors. Using the k values the 
transition lines defined by any of the observables can be given as 

Ufi) = T,[l - K ■ /T^] (4.9) 

with Tc once again being short for Tc(/i = 0). 

The results also suggest that the transition remains a weak crossover with essentially con- 
stant strength for small to moderate chemical potentials. Actually, there is a slight increase 
in the width of the transition determined from both quantities. This effect is, however, very 
weak: the width only changes by a few percent up to ~ T^. 
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Figure 4.7: The phase diagram of QCD for smah chemical potentials. The crossover transition 
between the 'cold' and 'hot' phases is represented by the colored area. The lower solid band shows 
the result for Tc(/i) defined through the chiral condensate and the upper one through the strange 
susceptibility. The width of the bands represent the statistical error of Tc(/x) for the given /i coming 
from the error of the curvature k for both observables. The dashed line is the freeze-out curve from 
heavy ion experiments [71]. Also indicated are with different symbols the individual measurements of 
the chemical freeze-out from RHIC, SPS (Super Proton Synchrotron) and AGS (Alternating Gradient 
Synchrotron), respectively. The center of mass energies ^/snn for each are shown in the legend. 

The validity range of our result is hard to estimate from the present study alone. A con- 
servative estimate for the limit where the result obtained through the Taylor-expansion is still 
reliable is /i ~ Tc i.e. where the expansion parameter exceeds unity. In the baryonic chemical 
potential this corresponds to about 500 MeV. Beyond this limit higher-order corrections are by 
all means expected to be important. Earlier experience with the exact reweighting method [50] 
also shows that the leading-order quadratic behavior of the Tc(/i) function dominates up to the 
above mentioned limit in the baryonic chemical potential. To investigate whether higher order 
terms may lead to a critical point one must carry out a similar analysis with full reweighting, 
beyond the reach of present computational resources. 

The final result is shown in figure 4.7. The crossover region's extent changes little as the 
chemical potential increases, and within it two definitions give different curves for Tc(/i). It is 
useful to compare the whole picture to the freeze-out curve [71] which summarizes experimental 
results on the {T, /i} points where hadronization of the quark-gluon plasma was observed. This 
curve is expected to lie in the interior of the crossover region, as is indicated by our results as 
well. 
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Chapter 5 



The equation of state with dynamical 
quarks 



As it was stressed throughout the introduction, the study of the finite temperature transition 
of QCD is relevant for both the case of the early Universe and that of heavy ion collisions. 
In both situations strongly interacting matter is thought to enter the high temperature QGP 
phase of quarks and gluons. While the nature and the temperature scale of the transition 
imply interesting cosmological consequences (see section 1.1), they also have a direct impact on 
contemporary high energy experiments in particle accelerators. In particular, the deconfined 
phase of QCD is claimed to have been produced in the ultrarelativistic heavy ion collision 
experiments at CERN SPS, RHIC at Brookhaven National Laboratory, ALICE at the LHC and 
the future FAIR at the GSI. The experimental results available so far show that the hot QCD 
matter produced experimentally exhibits robust collective fiow phenomena, which are well and 
consistently described by near-ideal relativistic hydrodynamics [72-74]. These hydrodynamic 
models need as an input an Equation of State (EoS) which relates the local thermodynamic 
quantities. Therefore, an accurate determination of the QCD EoS is an essential ingredient 
to understand the nature of matter created in heavy ion collisions, as well as to model the 
behavior of the QGP in the early universe. 

QCD thermodynamics and in particular the EoS have received increasing attention in the 
past years. These include among others calculations in the quenched approximation [75, 76], 
two-fiavor simulations [77, 78], studies with heavier-than-physical [79, 80], almost physical 
[69, 81-83] and physical quark masses [84]. However, results for an EoS involving physical 
masses together with a careful continuum limit are so far missing; their relevance for the 
physics of the hot deconfined matter is obvious. In section 5.3 I present results on various 
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thermodynamic observables for a system of A^g = 2 + 1 flavors of dynamical quarks. The 
findings have recently been published, see Ref. [12]. Two important aspects of this study 
are to be emphasized. First, the quark masses are set to their physical values, which is very 
important due to the strong dependence of thermodynamic quantities on these parameters. 
Second, a continuum limit estimate is given using lattices of temporal extent Nt = 6, 8, 10 
and a few checkpoints at A'^t = 12. Just as for the study of the phase diagram, once again the 
Symanzik improved gauge and the stout smeared staggered fermionic action is used. 

5.1 Determination of the pressure 

First it is instructive to review the standard integral method, which can be used to determine 
the equation of state on the lattice [85]. This technique is then further developed by the 
multidimensional spline method presented in section 5.1.2. 

5.1.1 The integral method 

The starting point is the (dimensionless) pressure as defined on the lattice 

p'^\f3,m,) = (iViiV3)-MogX(/3,m,) (5.1) 

as a function of the bare parameters, i.e. the inverse gauge coupling and the quark masses. 
This combination in itself is unfortunately not accessible using conventional algorithms, only 
its partial derivatives with respect to the bare parameters of the action are measurable. Using 
these partial derivatives the pressure can be rewritten as a multidimensional integral along a 
path in the space of bare parameters: 

p^-(/3,m,) -;.'-(/3o,m,o) = {N.N^r' ^"^'^ ( ^^^^ + E^^'^^^l ' (^.2) 

where the index '0' is used to denote the starting point of the integration. The value of the 
pressure at parameters '0' has to be handled with care: one either chooses the starting point 
so deeply in the strong coupling regime, so that the p''^*(/3o, "n^go) can already be neglected or 
it can be taken from model computations (see section 5.1.3). The derivatives in (5.2) are the 
gauge action density and the quark chiral condensate densities (as defined in (3.4)), measured 
in lattice units: 

(-^.) = (iv.ivi')-^^, = i^^^sr'^- (5-3) 
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As discussed in section 3.2.2, the renormalization of the pressure can be carried out by 
subtracting its value measured on a lattice with the same bare parameters, but at a different 
temperature value, ie. with a different temporal extent Nf^^. In the following for these 'cold' 
lattices we use Nf^^ = 3Nt, i.e. ^ = 3 in the notation of (3.29) and (3.32). The finiteness of 
N^^^ causes an error of the order which is smaller than the typical size of the statistical 
error of the results, see section 5.3. 

Using the subtracted observables 

one can express the renormalized pressure, divided by as: 

^ - ^ = iV.^ (d^i-s.r'' + Y: dm,{^P,^P,rA (5.6) 

where T and Tq are the temperature values corresponding to the lattice spacing at the bare 
parameters {/3,mg) and (/3o, "Ti^o)- Using (3.19) one can also relate the integrand to the trace 
anomaly: 

= (^dP{-s,r' + E dmSqi^.r''^ . (5.7) 

Let us make an observation here: due to the prefactor in the trace anomaly the subtracted 
observables decrease with the lattice spacing (at a fixed temperature this means an increasing 
Nf). While {—Sg)^^^ decreases as A^^"^, the chiral condensate behaves substantially different. 
Due to chiral symmetry^ ii'q'^q) is proportional to the bare quark mass, moreover it is also 
multiplied by another power of the bare quark mass (through the druq line-element) in the 
trace anomaly. Since the bare mass (in lattice units) also decreases with the lattice spacing, 
the subtracted condensate only decreases as Nf^. These scalings are directly translated into 
different precisions for the two terms in (5.7). The term with the gauge action density has 

times larger errors than the term with the chiral condensate, if they are evaluated with 
the same statistics. This is also reflected by the fact that in Sg all of the divergences of the 
free energy density remain, while in ipgipq the highest order term is 0(a~^). Therefore the 
subtracted observable (sg)^'^^ has larger cancellations as compared to {ipqipqY"^^ . 

The standard integral method proceeds as follows: first one calculates the trace anomaly for 
several temperatures along the lines of constant physics (LCP, see definition in section 2.4.1) 



^Staggered fermions have only a remnant chiral symmetry, but this does not affect the argument. 



59 



5. THE EQUATION OF STATE WITH DYNAMICAL QUARKS 



and then integrates it to get the pressure up to an integration constant (which can be either 
neglected or taken from a model calculation). This path was used in several lattice studies, 
e.g. [69, 82-84]. In the following a generalized version of the integral technique is presented. 



5.1.2 Multidimensional spline method 



In order to compare the pressure to results obtained with other fermionic formulations is is 
useful to study the light quark dependence of p. A major goal of the study presented in 
section 5.3 was therefore to determine the EoS for several (heavier than physical) different 
values of the light quark mass rriud = = f^d- To this end simulations were carried out 
covering an extended region in the {/3,ms,m„rf} parameter space. The strange mass was 
always set to its physical value m'^^^^{P) along the line of constant physics. On the other hand, 
rriud was set as rriudiP) = f^siP) / R with the quark mass ratio R ranging from 1 to Rp^'^'^ = 28.15. 
The value R = 1 corresponds to the three degenerate flavor case {Nq = 3), whereas Rp^^^ is the 
"real world" value of the quark mass ratio {Ng = 2 + 1). At these parameters we can consider 
the pressure as a function of only two variables: /3 and R. The respective derivatives can be 
measured on the lattice, they are given by the following combinations: 
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These derivatives have to be integrated to obtain 
the pressure. Since the integrand is in itself - by 
construction - a total derivative, the result has to 
be independent of the chosen integration path in the 
two-dimensional space of the parameters /3 and R. A 
few possible paths are illustrated in figure 5.1. The 
simulation points are denoted by crosses. A straight- 
forward way would be to perform a one-dimensional 
integral in /3 at a fixed value of the quark mass ra- 
tio R, this corresponds to path A. However one can 
also imagine zigzagging routes, like path B or C on 
the figure. Averaging several of such integrals, one 
can increase statistics and furthermore also provide 
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Figure 5.1: Illustration of possible inte- 
gration paths. 



an estimate of the systematic error related to the integration path. 
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Here I propose a generalized method that takes into account every possible integration path 
at the same time. The main idea is - instead of parameterizing the derivatives of the pressure 
and then integrate - to parameterize the pressure itself. A straightforward parameterization is 
one using the actual values pki of the pressure at some node points {/3k, Ri} that build up a two- 
dimensional spline function - i.e. that unambiguously determine the whole pressure surface as 
a function of /3 and R. The spline parameters pki can be set to minimize the deviation between 
the derivatives of this surface and the measured values D/s and -D^j. This minimum condition 
leads to a system of linear equations that can be solved for p^i. The solution for the spline 
parameters in turn determines the surface P{P,R). 

This method gives the pressure directly using the information contained in the derivatives 

and without the need to carry out an actual integration. Derivatives of the surface 
are then straightforward to compute, e.g. the trace anomaly can be obtained using the rela- 
tion (3.19). However, in order to ensure the smoothness of further derived quantities (e, s and 
Cs), a 4-degree spline fit to the trace anomaly was carried out to give the results in section 5.3. 
Using this method we are also capable of estimating the systematic error of the result (as op- 
posed to the conventional integral technique) by the variation of the nodepoints of the spline 
interpolation. Details of the spline fitting and the determination of the systematic errors can 
be found in appendix A and in [13]. 

Both the integral technique and the spline method can easily be generalized to the case of 
more quark flavors - of particular interest is the contribution of the charm quark. In order to 
obtain this contribution, the (3 derivative of the pressure in (5.8) has to be complemented by 
the charm condensate. With the charm mass set according to rric = Q ■ mg this derivative is 
modified as 

C/TTl — 

Dp^Dp + ^-Q- iiJci'cr'' (5.10) 
5.1.3 Adjusting the integration constant 

The multidimensional spline method developed in the previous section determines the pressure 
only up to an overall constant factor, which corresponds to an integration constant (this uncer- 
tainty is of course also present in the conventional integral method). This originates from the 
fact that a constant shift of the surface leaves its gradient unchanged - i.e. a shifted solution is 
an equivalently good fit to the measured derivatives. This integration constant unfortunately 
cannot be determined on the lattice since the derivatives Dj^ and Dji are measured at finite 
values of the temperature (where the actual value of the pressure is a priori unknown). In 
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order to obtain a final result for the pressure this constant shift needs to be estimated; here I 
present two possible approaches to do so. 

Set p to zero at small T in Nq = 3 theory 

A straightforward way to set the integration constant is to set the pressure to zero at small 
temperature, that is to say, somewhere deep in the hadronic phase where contributions to p are 
suppressed. This suppression is exponential in the hadron masses niH and the temperature: 
~ e~™^^. Therefore it is reasonable to set the zero at larger quark masses, where hadrons are 
also heavier and thus the contribution to the pressure is smaller. Taking into account these 
considerations we set the integration constant such that at the smallest temperature for three 
degenerate flavors (the largest rriud) the pressure is set to zero: 

p(T = 100 MeV, i? = 1) = (5.11) 

The error of this choice is hard to estimate from the lattice alone. According to the HRG 
model (see next subsection) however, the pressure at this point p^^*^(T)/T^ = 0.02 is much 
smaller than the typical statistical errors on the lattice. 

Set p according to Hadron Resonance Gas model 

Though it is possible to calculate the pressure solely using lattice methods, it is instructive 
to compare this result to model calculations, especially at low temperature, where lattice dis- 
cretization errors are expected to be large. Moreover, there is an additional effect that may 
influence the pressure in the hadronic phase. This is due to the fact that the staggered formu- 
lation does not preserve the flavor symmetry of continuum QCD, and as a consequence, the 
spectrum of low lying hadron states is distorted (see section 2.4.4). Recent analyses performed 
by various collaborations [33-35] have pointed out that this distortion can have a dramatic 
impact on the thermodynamic quantities. In order to quantify this effect, one can compare the 
low temperature behavior of the observables obtained on the lattice, to the predictions of the 
Hadron Resonance Gas (HRG) model. 

The HRG model has been widely used to study the hadronic phase of QCD in comparison 
with lattice data [86-88]. The model builds on the fact that the low temperature phase of QCD 
is dominated by pions. Goldstone's theorem implies weak interactions between pions at low 
energies, which allows to study them within chiral perturbation theory [89]. As the temperature 
T increases, heavier states become more relevant and need to be taken into account. In the 
HRG approximation [90], the microcanonical partition function of the interacting system can 
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be calculated in the thermodynamic limit V ^ oo, assuming that it is a gas of non-interacting 
free hadrons and resonances [91]. 

The pressure of the HRG model can therefore be written as the sum of independent con- 
tributions coming from non-interacting resonances 



^ J2 log Z'\T, V, f^x^,m,) + ^ log Z^(T, V, /ix-, m,) (5.12) 



j<4 yj'S — V- ) - ! r-^"! ■■"1/ ' VT^ 

iG mesons baryons 



where 



log Z'\T, V, fix^,mi) = I dke log(l - z.e-^»/^) 

Vd- 

log X^(T, V, fix'^,mi) = / dkk' log(l + z,e-^»/^) 

Jo 



(5.13) 



with energies Ei = ^/k"^ + mf, degeneracy factors di and fugacities 

Zi = exp(^^J2^>^^^ (5.14) 

In the above equation, X"- are all possible conserved charges, including the baryon number 
B, electric charge Q and strangeness S. The sums in (5.12) include all known baryons and 
mesons up to 2.5 GeV, as listed in the latest edition of the Particle Data Book [23]. The trace 
anomaly of the HRG model can also be calculated from (5.12) using the relation (3.19). 

As expected, for temperatures T < 60-100 MeV the contribution of the pions dominate the 
pressure. On the other hand, for larger temperatures the kaon contribution becomes sizable 
and slowly, heavier states also become relevant. The contributions coming from these sectors 
are illustrated by the colored bands in the left side of figure 5.2. One sees that at T ~ 100 
MeV (where a comparison to lattice results is desired) the pressure is almost completely due 
to the pionic contribution and resonances with large widths do not play a significant role. 

As already mentioned the staggered lattice discretization has a considerable impact on the 
hadron spectrum. In order to investigate these errors, we define a "lattice HRG" model, where 
in the hadron masses lattice discretization effects are taken into account. We consider only the 
taste splitting effects for the pions and the kaons. For example the contribution of the pions 
to the normalized pressure is now given by a sum, which runs over the different pion tastes: 

/iX",mQ,), (5.15) 

a=0 
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T[MeV] T[MeV] 

Figure 5.2: Contributions to the pressure at low temperatures according to the HRG model (left 
side). The normalized trace anomaly in the physical HRG model (solid red line) and in the HRG 
models with the lattice hadron spectrum (dashed blue lines) (right side). 

where tIq, are the multiplicities of the various taste sectors, and nia the masses of the pion 
tastes (see section 2.4.4). 

In the right side of figure 5.2 we plot the normalized trace anomaly of the HRG model 
with the physical and with the lattice distorted spectrum for the four different lattice temporal 
extensions used in our investigations. The difference between the physical and lattice curves 
is a first estimate of the lattice discretization errors arising from the taste violation. As 
the temperature decreases at a fixed Nt, the lattice spacing gets larger and so do the taste 
violation effects. The model calculation suggests, that the lattice results may have sizeable 
systematic errors in the low temperature region. Above T ~ 100 MeV, this error estimate for 
the interaction measure is smaller than the typical magnitude of other errors in lattice QCD 
calculations. 



5.2 Systematic effects 

As it was pointed out in chapter 2, there are two main sources of systematic error that one has 
to handle with care. One stems from the lattice artefacts and the other from finite size effects. 
In the next two subsections I will present results that support that neither of the above two 
systematic errors are present in our analysis of the EoS. 
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For high temperatures the thermodynamic quantities approach the corresponding values 
of the non-interacting massless relativistic gas, i.e. the Stefan-Boltzmann (SB) hmit. For the 
three flavor pressure the SB value is Psb/T'^ ~ 5.209. Using (3.21) one can also calculate the 
limit for the energy density ess = ^^Psb, the entropy density ssb = ^Vsb/T and the square of 
the speed of sound 55 = 1/3. 

In order to decrease lattice artefacts - besides using an improved action - we also carry out 
an improvement on the level of the observables. In lattice thermodynamics the cutoff effects in 
general depend on A^^^, however, for Nt ^ 00 they disappear. In the following we use a tree- level 
improvement for the pressure: we divide the lattice results by the appropriate improvement 
coefficients ^{Nt). At the tree-level these coefficients are proportional to the value of the 
pressure as measured in the noninteracting system on a lattice with finite temporal extent Nt: 

m) = (5.16) 

PSB 

For the action of our choice the n{Nt) factors are calculated to be: 



Nt = 6 


Nt = 8 


Nt = 10 


A^t = 12 


1.517 


1.283 


1.159 


1.099 



We use the same factor for the two different spatial volumes in our finite volume study (6 x 18^ 
and 6 x 36^). Using thermodynamic relations one can obtain these improvement coefficients 
for the energy density, trace anomaly and entropy, too. The speed of sound receives no im- 
provement factor at tree level. 

In the left side of figure 5.3 we illustrate at three temperature values (T = 132, 167 and 
223 MeV) the effectiveness of this improvement procedure. We show both the unimproved 
and the improved values of the normalized trace anomaly I/T^ for Nt = 6,8,10 and 12 as 
a function of 1/N^. The lines are linear continuum extrapolations^ using the three smallest 
lattice spacings. As visible in the figure, the continuum limit of both the unimproved and 
the improved observables agree well. Thus, our results confirm the expectations, that the 
lattice tree-level improvement effectively reduces cutoff effects. At all three temperatures the 
unimproved observables have larger cutoff effects (i.e. larger deviation from the continuum 
value) than the improved ones. 

Note that the improvement coefficients K,{Nt) are exact only at tree-level, thus in the in- 
finitely high temperature, non-interacting case. As we decrease the temperature, corrections 

^Note that the staggered fermionic action approaches the continuum action as ~ a^. This scahng motivates 
that the extrapolation is Unear in ^ a? . 
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to these improvement coefficients appear, which have the form 1 + b2{T)/N^ + .... Empirically 
one finds that the b2{T) coefficient, which describes the size of lattice artefacts of the tree- level 
improved quantities, is tiny not only at very high temperatures, but throughout the decon- 
fined phase. In particular, for the case of the three temperatures considered in the left panel 
of figure 5.3, the 62 (T) coefficients are found to differ from zero with less than one standard 
deviation. 
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Figure 5.3: The normalized trace anomaly at three different temperatures as a function of l/N^ 
Filled red symbols represent the results within the lattice tree-level improvement framework, blue 
opened symbols show the results without this improvement. The error of the continuum extrapolated 
value is about 0.4 for all three temperatures (left panel). The trace anomaly on lattices with different 
spatial volumes: red band for Ng/Nf = 3 and blue points for Ng/Nf = 6 (right panel). 



The second possible systematic error is that connected to the finite lattice size. In order 
to verify that there are no significant finite size effects present in the lattice data with the 
aspect ratio Ns/Nt = 3, we checked our Nt = 6 data against a set of high precision Ns/Nt = 6 
simulations. The latter lattice geometry corresponds to about 7 fm box size at the transition 
temperature. The right panel of figure 5.3 shows the comparison between the two volumes 
for the normalized trace anomaly I/T^. From this result we conclude that it is acceptable to 
perform the more expensive simulations throughout with Ng/Nf = 3. Let us note here, that 
the volume-independence in the transition region is an unambiguous evidence for the crossover 
type of the transition. 
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5.3 Results 

In this section I present results on the equation of state of Ng = 2 + 1 QCD. The pressure, the 
interaction measure, the energy and entropy density as well as the speed of sound are shown as 
a function of the temperature. The characteristic points of these observables are identified and 
a parametrization for the trace anomaly is presented. Afterwards I study the dependence of 
the pressure on the light quark mass rriud- Furthermore, the contribution of the charm quark 
to the pressure is discussed within the partially quenched framework. Finally a comparison to 
results obtained with other fermionic discretizations is carried out. 

The multidimensional spline method (see section 5.1.2 and appendix A) was used to de- 
termine the equation of state on A'^t = 6, 8 lattices for various different values of the quark 
masses. In the case of Aj = 10 we determined the equation of state exclusively with physical 
quark masses. In order to satisfy the normalization condition of (5.11), we made heavier mass 
simulations at T = 100 MeV up to the three flavor point. For A^ = 12 we made simulations 
at three temperature values, this allows us to calculate the trace anomaly only. 

The error bars on the figures are obtained by quadratically adding the statistical error and 
the systematic error of the spline method (see appendix A). The temperature values have an 
error at the 2% level arising from the scale setting procedure (see section 2.4.1). 



5.3.1 The Nq = 2 + 1 flavor equation of state 

In figure 5.4 the normalized trace anomaly is shown as a function of the temperature. As it 
can be seen in the figure, results show essentially no dependence on the lattice spacing, as all 
four datasets lie on top of each other. Only the coarsest Nt = 6 lattice shows some deviations 
around ~ 300 MeV. In the same figure, I also provide a zoom of the transition region; here 
the results from the HRG model are also plotted. A good agreement with the lattice results is 
found up to T ~ 140 MeV. 

One characteristic temperature of the crossover transition can be defined as the inflection 
point of the trace anomaly. This and other characteristic features of the trace anomaly are 
summarized in the following table: 

Inflection point of /(T) /T* I 152(4) MeV 
Maximum value of /(T) /T^ 4.1(1) 
T at the maximum of I{T)/T^ \ 191(5) MeV 

In figure 5.5 I show the main result of the present study: the pressure of A^^ = 2 + 1 QCD 
as a function of the temperature. Here results for three different lattice spacings are presented. 
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The Nt = 6 and Nt = 8 are in the temperature range from 100 up to 1000 MeV. The results 
on A'^t = 10 are in the range from 100 up to 365 MeV. The zero point of the pressure is set 
by (5.11), as supported by the discussion in subsection 5.1.3. Note that from this condition 
we get a nonzero pressure for the Ng = 2 -\- 1 system (i.e. R = 28.15, cf. (5.11)) even at the 
smallest temperature T = 100 MeV. 

The value of the pressure at 100 MeV is approximately two third of the value suggested by 
the HRG model. The origin of this difference cannot be clarified at the moment. One expects 
that the lattice artefacts are considerably larger at low temperatures, than what one estimates 
from the difference of A'^t = 6, 8 and 10 results. This is quite reasonable, since even our finest 
lattice at T = 100 MeV has ~ 0.2 fm lattice spacing, which is far from the regime, where 
lattice results starts to scale. On the other hand, this discrepancy might also point to the 
failure of the HRG model. In order to be on the safe side the size of this unexplained difference 
is considered as an estimate of our systematic uncertainty in the low temperature regime. 

In figure 5.6 I show the energy density (upper panel) and the entropy density (lower panel), 
while in the upper panel of figure 5.7 the square of the speed of sound as a function of the 
temperature. Finally, in the lower panel of figure 5.7 the speed of sound and the ratio p/e are 
presented as functions of the energy density. 

Some characteristic properties of the speed of sound and the ratio p/e are also of phe- 
nomenological interest. These are summarized in the following table: 
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T[MeV] 

Figure 5.5: The pressure normalized by as a function of the temperature on A'^t = 6, 8 and 10 
lattices. The Stefan-Boltzmann limit pssiT) ~ 5.209 • is indicated by an arrow. 



Minimum value of c^(T) 
T at the minimum of c^(r) 
e at the minimum of c^(T) 


0.133(5) 
145(5) MeV 
0.20(4) GeV/fm^ 


Minimum value of p/e 
T at the minimum of p/e 
e at the minimum of p/e 


0.145(4) 
159(5) MeV 
0.44(5) GeV/fm^ 



5.3.2 Continuum estimate and parametrization 

Based on the results for various lattice spacings we can also provide a continuum estimate 
^ for the quantities presented above. We take the average of the data at the smallest two 
lattice spacings and as an error we assign either the difference of the two or the statistical error 
depending on whichever is larger. As already mentioned, for low temperatures the lattice result 
for the pressure is significantly smaller than the prediction of the HRG model: at T = 100 
MeV the lattice result is p{T)/T^ = 0.16(4), whereas the model prediction is p{T)/T^ = 0.27. 
Therefore in the continuum estimate of the pressure we shift the central values of the lattice 
results up by the half of this difference (0.06) and this shift is then considered as a systematic 
error in the entire temperature range. 

The trace anomaly as a function of T plays a central role in hydrodynamic models and thus 

^For a rigorous continuum extrapolation one would need Nt = 12 for the entire temperature region. 
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Figure 5.6: The normalized energy density and entropy density as functions of the temperature on 
A'^t = 6, 8 and 10 lattices. The SB limits esB = '^PSB and ssb = ^Psb/T are indicated by the arrows. 



we find it useful to give a global parametrization that describes J/T^ in the whole temperature 
range considered. We considered the following fit function: 

liT) ( . „ . /,2^ (. , /o-[tanh(/i-t + /2) + l] ^ 

- exp(-/ii/t - /i2/t ) ■ Aio H T"^ — 75 h (5-17) 



\ 1 + gi-t + g2-V 

where the dimensionless t variable is defined as t = T/(200MeV). This function reproduces 
the continuum estimate for the normalized trace anomaly in the entire temperature range 
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Figure 5.7: The square of the speed of sound as a function of the temperature on A''^ = 6, 8 and 10 
lattices (upper panel) and and p/e as a function of the energy density on Nt = 8 lattices (lower 
panel). The SB limit of 1/3 for both quantities is indicated by the arrows. 

T = 100 . . . 1000 MeV. The {/o, /i, /2} parameters describe the steep rise of the trace anomaly 
in the transition region, whereas the {gi, (72} parameterize the decrease for higher temperatures. 
The parametrization also approximates the HRG model prediction for T < 100 MeV, this is 
described by the {/iq, hi, /12} parameters. The parameters can be found in table 5.1. 
For these temperatures the difference in the trace anomaly between the parametrization and 
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Figure 5.8: Continuum estimate for the trace anomaly normalized by together with the 
parametrization of (5.17) using the Ng = 2 + 1 parameters from table 5.1. 
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Table 5.1: Parameters of the function in (5.17) describing the trace anomaly in the A''^ = 2 + 1 and 
in the A''^ = 2 + 1 + 1 flavor cases. 

the HRG model is less than A(/(T)/T^) < 0.07. The parametrization together with our 
continuum estimate for the trace anomaly is shown in figure 5.8. From this parametrization 
the normalized pressure can be obtained by the definite integral of (3.20). The so obtained 
function goes through the points of the continuum estimate of the pressure for temperatures 
T = 100 . . . 1000 MeV and for T < 100 MeV the deviation from the HRG prediction is less 
than A{p{T)/T^) < 0.02. 

5.3.3 Light quark mass-dependence 

The multidimensional spline method gives the pressure (and the derived quantities) as a smooth 
function of both the temperature and the light quark mass. Thus it becomes possible to 
investigate the dependence of the EoS on niud- 

In figure 5.9, the trace anomaly for two different values of the light quark masses is plotted: 
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Figure 5.9: The normalized trace anomaly for two different values of the light quark masses on 
Nt = 8 lattices: the physical jn^^^ = ni^^^ and the three degenerate flavor = m^^-^^ case. 

for the physical case, where m^^ = Tir^^ and for the three degenerate flavor case, where 
f^ud = mf^^'^ (i.e. R = 28.15 and i? = 1 in the notation of section 5.1.2, respectively). This 
latter case corresponds to a pion mass of approximately ~ 720 MeV. The results are from 
our Nt = 8 lattices, this is the smallest lattice spacing, where we have the complete mass 
dependence of the equation of state. As it is expected, the peak position of the trace anomaly 
is shifted towards higher temperature values for larger quark masses. The position in the three 
degenerate flavor case is ~ 25% larger than at the physical point. The height also increases 
by about ~ 40%. When zooming into the transition region, we also show the comparison 
with the HRG model. For low temperatures one finds a reasonable agreement also in the 
heavy quark mass case. As it is expected, the dependence on the quark masses vanishes as 
one goes to higher temperatures. Therefore it is plausible to compare the result with massless 
perturbation theory. A good agreement can be observed with the highest order perturbative 
calculation without non-perturbative input {0{g^), see figure 5.9). 

5.3.4 Estimate for the Nq = 2 + 1 + 1 flavor equation of state 

While at low temperatures the equation of state only contains terms that originate from the 
light quarks, in high energy processes charm quarks can also be created from the vacuum, 
and they can also be present in the initial or final states. An interesting issue regarding the 
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EoS is whether the charm quark can give an important contribution to the QCD EoS, in the 
range of temperatures which are reached in heavy ion colhsions. It is often assumed that it 
can be neglected, the charm mass being too heavy to play any role at T ~ 2Tc. However, 
perturbative QCD predicts that its contribution to thermodynamic observables is relevant at 
surprisingly low temperatures, down to T ^ 350 MeV [92]. Recent exploratory lattice studies 
have confirmed these expectations [93, 94], indicating a non-negligible contribution of the 
charm quark to thermodynamics already at 1.5Tc. These results have been obtained on rather 
coarse lattices {Nt = 4, 6) and with the charm quark treated in the quenched approximation. 
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Figure 5.10: Contribution of the charm quark to the pressure for two different lattice spacings 
(upper panel). The pressure normalized by T'^ for Ng = 2 + 1 + 1 and Ng = 2 + 1 flavors on A^^ = 8 
lattices (lower panel). The corresponding Stefan-Boltzmann limits are indicated by arrows. The 
charm to strange quark mass ratio is Q = 11.85 on this plot. 



We estimated the contribution of the charm quark to the pressure on our Nt = 6 and 8 
lattices at the partially quenched level (i.e. we used the same configurations as for the Ng = 2+1 
case). The measurements were carried out for several values of the charm to strange quark mass 
ratio Q = mc/mg. According to a recent high-precision lattice calculation [95] the physical 
value of this ratio is Q^^^^ = 11.85(16). For this central value we show the contribution Pcharm 
as a function of the temperature on the upper panel of figure 5.10. We find that Pcharm is 
non-zero already at temperatures T ~ 200 MeV. The total A'g = 2 + 1 + 1 pressure is compared 
to the already presented A^^^ = 2 + 1 pressure on the lower panel of figure 5.10. Using the 
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parametrization in (5.17) we also fit the Nq = 2 + 1 + 1 data for the trace anomaly for the 
Nt = 8 data. The resulting fit parameters can be found in table 5.1. 

It is important to stress here that the estimate for the charm contribution presented here 
suffers from two uncertainties. First, we neglected the back-reaction of the charm quarks on 
the gauge field as the charm condensate was measured at the partially quenched level. Second, 
due to the large mass of the charm, large lattice artefacts can also be expected. In order to 
provide a more precise calculation of the A'^ = 2 + 1 + 1 flavor equation of state, one has to relax 
the uncontrolled partial quenched approximation and introduce the charm quark dynamically. 
Moreover, the lattice spacing has to be further reduced to ensure a good resolution of the 
charmed excitations. The reduction of lattice artefacts can also be done by using improved 
fermion actions, therefore these are also increasingly important when studying heavy quark 
degrees of freedom.^ 

5.3.5 Comparison with different fermion discretizations 

In the past years, the equation of state of QCD was studied with various different fermion 
discretizations. Among these are improved staggered actions like the p4, the asqtad and the 
hisq action (see section 2.4.5). The disagreement in thermodynamic quantities between the 
previous two formulations and the stout smeared action used in the present thesis has a broad 
literature [24, 33, 35]. The main difference can be described by a ~ 20-30 MeV shift in the 
temperature. In particular, the transition temperatures obtained from various thermodynamic 
observables like the renormalized chiral condensate and the quark number susceptibilities differ 
significantly. 

It is interesting to look for this discrepancy in the equation of state as well. In 5.11 I 
show our results for the trace anomaly and a comparison with results obtained the p4, the 
asqtad and the hisq action. The data for the latter three are taken from [69], [83] and [96], 
respectively. As it can be clearly seen, the upward going branch and the peak position are 
located at ~ 20 MeV higher temperatures for the p4 and the asqtad actions as compared to 
the stout results. Furthermore, the height of the peak is also about 50% larger for the former 
two. One can also note that a weaker disagreement is observed between the stout and the 
hisq results. The origin of this discrepancy is most probably connected to lattice discretization 
errors and artefacts resulting from the taste splitting of the staggered tastes. Accordingly, 

^We remark that preliminary results regarding the dynamical contribution of the charm quark indicate that 
the above results may overestimate the iVg = 2 + 1 + 1 pressure. 
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Figure 5.11: The normalized trace anomaly is compared to recent results obtained with different 
fermionic actions [69, 83, 96]. 

results with both actions, as measured with smaller lattice spacings and smaller splitting seem 
to approach each other so that this discrepancy may be resolved in the near future. 
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Chapter 6 



The equation of state at high 
temperatures 

Although results from RHIC support the idea that the high temperature QGP behaves like a 
strongly coupled liquid, there are observed phenomena (like jet quenching or elliptic flow) that 
can be described well also by perturbative methods (i.e. with weak couplings). The question, 
whether the deconfined phase of quarks and gluons is more like a weakly coupled plasma or 
a strongly coupled liquid is thus of high relevance. A possible way to test the validity of 
perturbative approaches is to make a comparison with non-perturbative results obtained from 
the lattice. However, while standard perturbation theory converges only at extremely high 
temperatures, available lattice results end at around (5 — 10) ■ Tc. Including higher order terms 
in the expansion and applying resummation techniques can certainly push the lower limit of 
the validity of perturbation theory to smaller temperatures. Nevertheless, in order to carry out 
a proper comparison to perturbation theory, lattice measurements also have to be extended 
to higher temperatures. In the pure gluonic sector this has become recently possible. In this 
chapter I show results obtained in the SU{3) theory with the Symanzik improved gauge action 
up to previously unreached temperatures and also make a comparison to various perturbative 
approaches. 

6.1 Perturbative methods 

Perturbation theory has now been long used to calculate thermodynamic quantities in the 
weakly interacting (high temperature) limit. Most importantly, the pressure has been studied 
and recently determined up to order 0{g^\og{g)) [97]. Apparently, perturbative precision 
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cannot be pushed further to higher loops due the appearance of serious infrared divergences. 
These divergences are in general characteristic to finite temperature field theories. In particular, 
these problems prevent one to calculate the pressure in the order 0{g^), namely there is an 
unknown coefficient in this contribution. Unfortunately, even computed to this high order, 
the expansion only converges at extremely high temperatures, since at low temperatures the 
various terms oscillate strongly and thus can produce large cancellations. 

The disagreement between lattice data and standard perturbation theory has also been 
studied for other quantities like the expectation value (^^) at large /3 i.e. small lattice spac- 
ings [98]. The deviation from non-perturbative results apparently stems from the fact that 
the bare coupling constant is used as an expansion parameter in the perturbative series. One 
could very well define the perturbative expansion also in terms of a renormalized coupling con- 
stant Qr, which can be defined using some physical observable. Therefore, Qr will be a running 
coupling QrifJ'), which depends on the relevant energy scale i.e. the renormalization scale /i^. 

From a different point of view, the reason for the poor convergence of perturbative series 
can also be due to various plasma effects that influence the high temperature description 
of the system. These effects can be included by means of a systematic resummation of the 
perturbative series, like in e.g. the Hard Thermal Loop (HTL) approach [99]. Whether such 
resummed expansions do reproduce lattice results at reasonably small values of the temperature 
is an open question and is one of the main aims of our study. 

Moreover, there are also indications that there is a non-perturbative contribution to the 
equation of state that can never show up in perturbation theory [100]. For dimensional reasons, 
any flnite order perturbative formula will only give logarithmic corrections to the p{T) ~ 
Stefan-Boltzmann law. Instead of such logarithmic corrections, lattice data suggests that there 
is a two-dimensional "condensate", proportional to T^, which gives a dominant contribution 
for temperatures up to ~ 4Tc. This can be best observed in the trace anomaly I = e — 3p 
(see deflnition in section 3.1.4). Speciflcally, it is instructive to study the combination J/T^ ■ 
(T/Tc)^, which is shown in flgure 6.1. Our results with the Symanzik improved gauge action 
for various lattice spacings are compared to results obtained with the Wilson gauge action [75], 
the 5-loop perturbative expansion [97] and the HTL NNLO scheme. While for the former the 
renormalization scale fi = 27tT is used (black dashed line in the flgure), for the latter a range 
of f^uTL = ttT . . . AttT is considered (gray band) . 

^The former identification of n with the chemical potential is abandoned from now on. 
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Apparently, the combination J/T^ as mea- 
sured on the lattice is approximately constant 
in the temperature range Tc < T < 5Tc 
(there are however discrepancies between the 
Symanzik and Wilson results, see discussion 
later). While up to 5Tc lattice results seem 
completely incompatible with the perturba- 
tive predictions, at larger temperatures our 
results also account for the T^'-like steep rise 
in /(T) indicating a qualitative agreement 
with perturbative methods. In view of the 
fact that the trace anomaly measures the de- 
viation from the ideal gas (see discussion in 
section 3.1.4), this suggests that besides the 
ideal (perturbative) contribution ~ T^, / also 
contains a non-ideal (non-perturbative) term 
~ T^. Thus we separate the trace anomaly 
into two parts: 
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Figure 6.1: Our results for the normalized trace 
anomaly multiplied by T^/T^ for Nt = 5,6,7 and 
8 (red, green and blue dots, respectively). Also 
plotted are lattice results of [75], perturbation 
theory [97] and the HTL approach [99]. 
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The pressure can be obtained from I/T^ as a definite integral (see (3.20)). At extremely high 
temperatures its value is given by the Stefan-Boltzmann limit psB = 87r^/45T^. Integrating 
down from this point one obtains 
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The results for the trace anomaly in the high-temperature region allow for a fitting of the 
HTL renormalization scale /irtl and the unknown coefficient {qc in the notation of [97]) of 
the 0{g^) order contribution of perturbation theory. While Qc has already been calculated by 
means of a fit to the lattice data of [75], here we are able to repeat this fitting procedure at 
a much higher temperature, where the sixth order can be shown to be a minor correction to 
the fifth order. Once the optimal coefficient of the term is known, the non-perturbative 
contribution can also be quantified through a fit to some specified parameters of the function 
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We remark that the observed non-perturbative pattern may be explained within a fuzzy bag 
model [100], in terms of a dimension-2 gluon condensate [101, 102], in a system of transversely 
polarized quasi-particles [103] or within the gauge/string duality [104]. Here we do not go into 
the viability of such models and only identify it as the dominant non-perturbative contribution. 

6.2 Finite volume effects 

Existing lattice results for the pressure end at around ST^. These include results in the pure 
gauge sector with the Wilson plaquette action [75, 105] and also with various improved actions 
like the Symanzik action [106, 107], renormalization group-improved actions [76] or fixed-point 
actions [108]. The effect of changing the number of colors [109] was also studied. Results for 
the pressure of full QCD [79, 81, 84] are also present only up to (5 — 10) ■ Tc, like for the study 
presented in the previous chapter [12]. 

There are two main reasons for the absence of high temperature results: first, at increas- 
ingly high temperatures the signal/noise ratio in the trace anomaly decreases significantly. 
Consequently, it becomes more and more difficult to detect a nonvanishing value for J/T^. 
Note that this small signal is just the information necessary to fit the unknown perturbative 
parameters, as discussed in the previous section. (This is the reason that the primary observ- 
able in our study is the trace anomaly and not the pressure, as in chapter 5). Second, since 
the lattice spacing varies with the temperature as a = {NtT)~^, in order to have a constant 
physical lattice size , the number of lattice points in the spatial directions A^^ in principle has 
to increase like T. While the former problem can be avoided by accumulating larger statistics, 
the latter obstacle is more of a matter of principle. 

There is, however, a perturbative argument that may help one to improve on the situation. 
In the high temperature region the relevant correlation length is the scale that corresponds 
to the electric and magnetic screening masses. These are given by ~ gT and ~ g^T, 
respectively, and have been calculated on the lattice [110]. Thus, according to the magnetic 
mass (which is the one that gives larger correlations), in the high temperature region one has 

N,a-mrn>l ^ NjNt>T/mrn^ P (6.3) 

which implies that the aspect ratio has to be increased only proportionally to /3 = Q/g"^. 
However it is important not to overlook the loophole in this argument. On the lattice at a 
given temperature, one can never know whether the perturbative region is already reached. 
Justifying the neglection of finite size effects based on a perturbative reasoning, therefore, is 



80 



6.2 Finite volume effects 



not really convincing. In order to establish the range of validity of the perturbative approach 
itself, it is indispensable to simulate the non-perturbative scale ~ Tc, too. In this sense the 
aspect ratio sets the maximum temperature where the system is still non-perturbative as 
T < T, ■ NjNt. 

Furthermore, it is well known that the finiteness of the spatial volume can also distort the 
p = —f relation, as infrared corrections may appear. The impact of these corrections on the 
equation of state as a function of the aspect ratio has been analytically calculated recently [111]. 
According to this result, finite volume effects in the free energy density only drop below one 
percent if Ng/Nt > 8 holds for the aspect ratio. 

Keeping in mind the above considerations we perform three sets of simulations. First, we 
calculate the trace anomaly in the temperature range of T/Tc = 0.7 ... 15 (on 80^ x 5, 96'^ x 6 
and 112^ X 7 lattices) and also extract a continuum limit from these results. These lattices with 
aspect ratio Ng/Nf = 16 accommodate the non-perturbative scale up to approximately 
16Tc. We also support the continuum extrapolation with an additional Nt = 8 set of lattices 
(64^ X 8) below 8Tc. This combined extrapolation is described in the beginning of section 6.4. 
From the trace anomaly various other thermodynamic functions can be determined according 
to the relations of section 3.1.4. 

In the next step we study the finite volume scaling of the trace anomaly on a non-continuum 
data set at Nt = 5. We presents results using lattices of aspect ratio Ng/Nt = 4,6,8, 16 and 
24. The latter 120'^ x 5 lattice accommodates the scale up to 24Tc. Using these results we 
test finite size effects in the whole temperature region, see subsection 6.4.2. 

Since due to computational limitations we cannot increase the aspect ratio further, in our 
third set of simulations we calculate the continuum equation of state in a small box (on 40^ x 5, 
48^ X 6 and 64^ x 8 lattices) up to lOOOTc. In subsections 6.4.3 and 6.4.4 we use this data 
set to find the optimal free parameters of existing perturbative calculations, i.e. the already 
mentioned Qc parameter of 0{g^) perturbation theory and the renormalization scale /xhtl of 
the HTL scheme. Using these small volume results we observe a good agreement with the 
newly fitted perturbative formulae, indicating that this approach successfully connects the 
low temperature non-perturbative region with the high temperature perturbative realm. The 
precision of our data points exceeds any previous calculation by about an order of magnitude. 

In order for the large lattices to fit in the memory of our computer system, the renor- 
malization of the trace anomaly was done via half-temperature subtraction, as explained in 
section 3.2.2. Specifically, we calculate 

/(T) _ fI{T) 1 /(T/2) \ 1 /(T/2)-/(0) 
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Due to the fact that the trace anomaly very strongly depends on the volume around the 
critical temperature, in order to account for this dependence we measured the second term of 
the above expression on lattices with half the spatial extension. Thus both terms are measured 
with the same physical volume which ensures that the sum is smooth around T = 2n -Tc for 
each n G IN"*" (otherwise the difference between the volumes shows up as small cusps at these 
temperatures). For the lattices with half the spatial size (i.e. the second term in (6.4)) the 
subtraction is carried out in the standard way, i.e. at T = 0. The continuum limit from this 
combined technique is equal to what one finds using the standard scheme. 



6.3 Scale setting 

Besides the proper treatment of finite volume effects another challenging issue was the accu- 
rate determination of non-perturbative beta function corresponding to the Symanzik improved 
action. Instead of setting the scale using a or tq, it was advantageous to define the lattice 
spacing in terms of the transition temperature, as it was argued for in section 2.4.2. To this 
end we determined the critical couplings (3c up to Nt = 20 from the peak of the Polyakov loop 
susceptibility. We define the improved coupling in the "E" scheme [112], generalized for the 
Symanzik gauge action as 



gUg' 



2x1 



di di (6.5) 

Sl'~'\g') = d,g' + d^g' + . . . 

with coefficients cq and ci as introduced in (2.33). 

The expansion coefficients di we calculated explicitly using measurements of Pixi and P2XI 
at high /3, and checked against existing results in the literature [113]. We then matched the 



measurements of /3c to the universal two- loop running and converted to the MS scheme [114]. 
In figure 6.2 we plot as a function of N^'^ the lambda parameter in units of the transition 
temperature. We extrapolate to the continuum limit using a linear and a logarithmic fit, 
and define the systematic error of the result as half of the difference between the two. As a 
comparison we also show results for Tc/Ajj-^ in the lattice scheme, where the convergence is 
expected to be worse. 

Our final result we quote as 

' = 1.26(7) (6.6) 



82 



6.4 Results 



The error is overwhelmingly systematic 
and reflects the sensitivity to various contin- 
uum extrapolations. This is consistent with 
the combination of previous determinations: 
the Lambda parameter Ajjg = 0.614(2) (5)ro 
of [115] can be translated to a units using 
aro = 1.192(10) (based on [116]) and then 
used with Tja = 0.629(3) of [75]. Through 
our direct result one can easily translate the 
scale setting of the perturbative expressions 
to the lattice language. 
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Figure 6.2: Setting the scale by critical couplings 
in the E scheme. 



6.4 Results 

First we reproduce the results of [75] in the transition region. In figure 6.3 these results are 
compared to the trace anomaly measured on our first set of simulations, i.e. on large lattices 
{Ns/Nt = 16) with Nt = 5, 6, 7, supplemented by Nt = 8, with Ng/Nt = 8. From these four sets 
of results we perform a continuum extrapolation via a combined spline fitting method. The 
datasets for different lattice spacings are fitted together by an iVj-dependent spline function^. 
This "multi-spline" function - defined upon a set of nodepoints Pk with k = 1 ... K - is 
parameterized by two values at each nodepoint, written in the form + b^Nf'^. We fit these 
2K parameters to the measurements: the minimum condition for leads to a set of linear 
equations, which can be solved for the parameters. 

As a result we have a smooth function interpolating our data for each Nt (colored lines in 
the figure), together with a smooth, continuum extrapolated curve (yellow band in the figure). 
The statistical error of the fit is determined by a jackknife analysis, while the systematic error of 
the continuum result is estimated by the difference between the extrapolation from Nt = 5, 6, 7 
and Nt = 5, 6, 7, 8. Note that this method ensures an exact N^"^ scaling at each temperature. 
As visible in the figure, data points for various lattice spacings are on top of each other, with 
the exception of the transition region. This region is zoomed into in the inlay of the figure, 
showing that our data indeed exhibits the expected N^"^ scaling. 

^The A^f-dependcncc is of the form -/Vj~^, as motivated by the scahng properties of the Symanzik improved 
action. 
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Figure 6.3: The trace anomaly for various lattice spacings in the transition region. The result of a 
combined spline fit for each lattice spacing, together with the continuum extrapolation is shown by 
the colored lines and the yellow band, respectively. For comparison results with the standard Wilson 
action [75] are also shown by the dashed-dotted line. 

As figure 6.3 shows there is an apparent discrepancy between our continuum result and 
that of [75], particularly around Tc. This might be attributed to the finite volume effects as 
well as to the difference in the scale setting procedure. 

6.4.1 Comparison to the glueball gas model 

In order to explore the thermodynamics of the confined phase, next we zoom into the low 
temperature region T < Tc. In this region one can also calculate the trace anomaly within the 
glueball resonance model. In figure 6.4 we plot our results together with the contribution of the 
first twelve glueballs of [117]. There is an apparent deficit of the model prediction as compared 
to the lattice results. It has been suggested that this difference might be reduced if one allows 
for a temperature dependence of the glueballs [118], which has already been determined for 
the 0"'""'" and 2+"*" states [119]. This point was raised in connection to the lattice results of [109] 
below Tc. Our data is in fairly good agreement with this scenario, however further studies are 
necessary to understand the question in more detail. 

The glueball resonance model can also be used to set the integration constant in the pres- 
sure. We appoint as a reference temperature T^gf = 0.75Tc, and shift the lattice results up to 
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Figure 6.4: The trace anomaly in the confined phase measured on lattices with various lattice 
spacings and the continuum extrapolation (yellow band). The dashed line corresponds to the glueball 
resonance model, estimated from the twelve lightest glueballs. 

the resonance gas prediction at this temperature. This shift causes a ~ 0.1% change in our 
high temperature results for the pressure and is smaller than our statistical errors here. This 
reference point is used also for the pressure as measured in our third set of simulations (i.e. on 
lattices with N^/Nt = 8). 

6.4.2 Volume dependence of the results 

The effect of the non-perturbative ~ contribution to the trace anomaly reduces at increasing 
temperatures. Moreover, the presence of this contribution becomes unnoticeable at sufficiently 
high T, regardless of whether or not the lattice size accommodates the inverse scale. One 
way to discuss the relevance of this non-perturbative scale is to compare the trace anomaly 
at various spatial volumes. This comparison is shown in figure 6.5 for our Nt = 5 lattices. 
The standard aspect ratio Ng/Nt = 4 gives somewhat smaller values for I/T^, but beyond 
Ng/Nt = 6 we do not see any difference in the results above the transition region. 
Based on our observations so far we conclude that: 

1. Our large volume results for the trace anomaly agree well with the perturbative formulas 



above lOTc. 
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Figure 6.5: Volume dependence of the trace anomaly on our Nt = 5 lattices. 

2. Within statistical accuracy we see no discrepancy between results from different volumes 
provided that N,/Nt > 6. 

3. The significance of the non-perturbative contribution that is dominant in the transition 
region reduces as ~ 

4. All of our volumes accommodate the perturbatively relevant scales, in particular the 
magnetic screening length ~ [g'^T)~^. 

These considerations suggest that - even if the lattice volumes are ever shrinking as the tem- 
perature is increased - our results are able to describe the physical trace anomaly within the 
error bars shown. Thus it is reasonable to conjecture that our Ng/Nt = 8 dataset reliably 
connects the transition region with the perturbative regime. 

6.4.3 Fitting improved perturbation theory 

Regardless of whether the conjecture of the last subsection is valid or not, we can make use of 
our small volume simulations at high temperature to compare to perturbative expansions, in 
particular, to extract some unknown coefficients of these formulas. We perform the continuum 
extrapolation in the same manner as for the large volume data, see section 6.4, using the 
A'^t = 5, 6 and 8 lattices. The systematic error we define in this case as the difference between 
the continuum extrapolated and the Nt = 8 results. 
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First we compare our results to 0{g^) improved perturbation theory [97]. We perform a fit 
of the subtracted trace anomaly 

ipert{T, Qc, /X) - Zpert(T/2, Qc, /i) (6.7) 

to the unknown coefficient Qc of the term with a fixed renormalization scale of /i = 27iT. If 
we also allow for a variation of the scale, we find ^/2ttT to be consistent with 1 within errors. 
These fits are carried out for our results between lOTc < T < lOOOTc and the systematic error 
is estimated by varying the endpoints of the fit interval. Beyond this we also consider as a 
source of systematic error the uncertainty of our lattice scale setting (6.6). We quote as our 
final result for these parameters 

gc = -3526(4) (55) (30) (6.8) 

with the numbers in the parentheses are from left to right the statistical error, the error 
coming from the lattice scale and that from the variation of the fit interval. A good fit quality 
is indicated as x^/dof = 0.7. The fitted function is shown in figure 6.6. 




Figure 6.6: The continuum limit obtained from the lattice results (red band), compared to fitted 
perturbation theory. We fit the coefficient (gray dashed-dotted line), and the non-perturbative 
contribution in two different forms (black dashed line and green dotted line). 

While the ~ T^^ behaviour of the trace anomaly in the low-temperature region has been 
seen in many studies, e.g. [75, 109], its relative weight in the total observable has not yet been 
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quantified. Therefore we also consider it useful to estimate the non-perturbative contribution 
to the trace anomaly, which we assume to be of the form inp(T) = a^p + 6np log(T/Tc), i.e. we 
propose the following fit function: 

V.(r) + "'-^'''MTm (6.9) 

First we perform the fit for with b^p = kept fixed, then we carry out the fit for both 
non-perturbative coefficients. The fit interval is chosen to be < T < GOT^. We find that 
the constant approximation is not able to resolve the trace anomaly in the low temperature 
region as x^/dof ~ 150. The logarithmic correction improves much on the situation and we 
get x^/dof = 4.7. Moreover, the parameters are rather sensitive to the variation of the lower 
endpoint of the fit interval. Nevertheless, since there is no a priori constraint on the form of the 
fit function (6.9), we accept this as a first approximation to the non-perturbative contribution. 
We obtain the following coefficients: 

const.: flnp = 0.879(2) (40) &np = 

const. +log.: a^p = 1.371(1)(50) b^p = -0.618(2)(4) 

with the errors coming from the statistics and the lattice scale, respectively. We also plot both 
the constant and the logarithmic fit in figure 6.6. 

Using (6.2) the fitted perturbative formulae for the pressure are also straightforward to 
write down. In figure 6.7 we compare our small volume results to the so obtained predictions. 
Similar comparisons can be made for the case of the energy density and the entropy density 
also, where we find qualitatively the same behavior as for the pressure, see figures 6.8 and 6.9. 



6.4.4 Fitting HTL perturbation theory 

Next we discuss the region of validity of the HTL resummed perturbation theory. In particular, 
we compare once again our Ng/Nt = 8 lattice results to the NNLO expansion of the HTL 
scheme [99]. We consider the renormalization scale ;Uhtl as a free parameter of this expansion, 
and perform a fit to this parameter, i.e. our fit function to the subtracted trace anomaly is 

ipert(T, fiuTh) - ipcvt{T/2, /xhtl) (6.10) 

The fit is carried out for T > lOOT^, and the endpoint is varied to obtain the systematic error 
coming from the fitting procedure. The sum of deviations for this fit is x^/dof = 0.6, indicating 
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Figure 6.7: The normalized pressure as measured on our small volume boxes. A comparison is 
shown to various fitted perturbative functions. 
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Figure 6.8: The normalized energy density as measured on our small volume boxes. 



a nice agreement between lattice results and the perturbative expansion. Our result for the 
renormalization scale is (in the same notation for the errors as before) 

/^HTL 



27rT 



1.75(2)(6)(50) 



(6.11) 



The fitted formula for the trace anomaly is shown in figure 6.10. 
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Figure 6.9: The normalized entropy density as measured on our small volume boxes. 




1 10 1 00 



Figure 6.10: The trace anomaly measured on two different spatial volumes (red and blue 
bands), compared to the NLO and NNLO HTL expansion with varied renormalization scale 0.5 < 
/^HTL/2'/rr < 2 (green and gray shaded regions). The dashed-dotted line represents the expansion 
with the fitted scale of (6.11). 
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Chapter 7 



Summary 



In this thesis I studied the finite temperature transition between confined and deconfined matter 
at zero and nonzero quark densities. The findings are relevant for the understanding of the 
evolution of the early Universe and contemporary and upcoming heavy ion experiments. Our 
results were obtained using large-scale simulations using lattices with various physical extents 
and lattice spacings, at physical values of the parameters of the theory. In all three projects we 
thoroughly checked finite size effects and discretization effects, thereby controlling the two most 
important possible sources of systematic error. In order to ensure that the extrapolations lead 
to the correct continuum limit, all the necessary additive and multiplicative renormalizations 
were applied to the studied observables. 

In chapter 4 the study of the phase diagram of QCD at small chemical potentials was 
presented. We determined the curvature of the transition line which separates the chirally 
symmetric and asymmetric phases using two observables that are sensitive to this transition. 
Our results also contain information about the relative change in the width of the transition. 
This suggests that the strength of the transition decreases very slightly, therefore the weak 
crossover that is present at zero density persists up to quark chemical potentials fiu,d ^ Tc. 

In chapter 5 I showed results regarding the equation of state in full QCD. For this project 
I developed a new method to determine the pressure from the raw lattice data. This approach 
improves on the conventionally applied integral method such that it better estimates the sys- 
tematic error of the result, and is also capable of quantifying the quark mass dependence of 
the pressure. From the pressure we derived various other thermodynamic observables, which 
can be directly used to set the equation of state in e.g. hydrodynamic models. We estimate 
the effect of a charm quark in the pressure and show that it is significant already at T ~ 2Tc. 
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Furthermore, we observe that the pressure, the energy density and the entropy density are 
roughly 15-20% under their Stefan-Boltzmann hmits at our highest temperature ~ QT^. 

A these values of the temperature a reliable comparison with perturbation theory cannot 
be carried out since such perturbative expansions are only convergent at much larger temper- 
atures. In chapter 6 I report results of the equation of state in pure gauge theory at extremely 
high temperatures ~ lOOOTc. We perform simulations using lattices of size up to N^/Nt = 24 
and gather statistics that exceed results in the literature by at least an order of magnitude to 
determine the high temperature signal in the interaction measure. At large T we successfully 
match our results with perturbation theory and we are able to fit the unknown 0{g^) coeffi- 
cient of the expansion. We also quantify the non-perturbative contribution to the interaction 
measure, which is proportional to T^. 
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Appendix A 

Multidimensional spline integration 



In this appendix I describe the method, which was used to obtain the pressure as a function 
of the parameters /3 and R using its measured derivatives (5.8) and (5.9). This method is not 
exclusively applicable to the case of the pressure in QCD; on the contrary, it can be used in a 
generalized context, e.g. for statistical physical problems, where using standard Monte-Carlo 
methods one is not able to measure the (logarithm of the) partition function log Z itself, but 
only its derivatives with respect to the parameters of the theory. 

Therefore in the following discussion I introduce the algorithm for the general case, for 
reconstructing a multidimensional surface, using the gradient of the surface measured at some 
values of the coordinates. The algorithm includes an introduction of a set of nodepoints, 
upon which the multidimensional spline is determined. This determination is linear and thus 
straightforward to compute. By a systematic variation of the number and position of the 
nodepoints the method is then adapted to the particular problem. Unlike a multidimensional 
integration along some path, the present method results in a continuous, smooth surface, 
furthermore, it also applies to input data that are non-equidistant and not aligned on a rect- 
angular grid. Function values, first and second derivatives and integrals of the surface are easy 
to calculate. Moreover, the proper estimation of the statistical and systematic errors is also 
incorporated in the method. 

First I briefly summarize how a spline function can be introduced in an arbitrary number 
of dimensions D. Afterwards I show how the fit to the measurements is carried out. Since in 
practice D > 2 is seldom necessary, in order not to complicate the notation, the method is 
discussed in detail for the case of two dimensions. Nevertheless, D > 2 is also straightforward 
to implement. Finally I present the algorithm for the systematic placement of the nodepoints. 
The present method was recently published [13]. 
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A.l Spline definition in arbitrary dimensions 



A detailed introduction to the theory of sphnes can be found in [120-122]. For our purpose it 
suffices to know that in e.g. two dimensions a cubic sphne is defined upon a grid {xk, yi} with 
0<A;<-R',0</<L. The sphne surface is unambiguously determined by the values that 
it takes at the nodepoints fk^i = S{xk,yi) (and the boundary conditions, which we take to be 
"natural"). A grid square [xk,Xk+i] x [yi^yi^i] will be shortly referred to as {k,l}. The spline 
function itself is compactly written as^ 



3 3 

S{x,y) = ^ ^ C'rt/^ (fc) 

i=0 j=0 



if {x,y) G {k,l} 



(A.l) 



where Cfj are the spline coefficients and and 
the dimensionless coordinates: 



tik) = 



with the widths 



X Xk 



w 



y-yi 

(y) 
w)' 



(A2) 



(y) 

W' = yi+i 



yi 



(A.3) 



L-2 
2 

1 



1=0 
k=0 



K-2 K-1 



The linear equations that determine the spline param- 

eters fk,i from the coefficients for the x and y directions Figure A.l: Grid points in two dimen- 
(summarized in matrix form as X and F, respectively) sions. 

are independent from each other. Inverting these equations one obtains for the coefficients 



K-l L-1 



c. 



k,l 



l\ni £ 
) Ak+iJ ni,n2 



(A.4) 



ni=0 n2=0 



Here the matrix X {Y) only depends on the number K (L) of grid points in the x {y) direc- 

Similarly, in arbitrary dimensions D the different directions 



tion and the widths w 



{^) („,Xy) 



w 



decouple and thus (A.4) is easily generalized to the case of more dimensions. In the following 
we introduce the fitting method in two dimensions; the generalization for higher dimensions is 
also straightforward to implement. 



^Note that in this formulation the cubic sphne contains terms hke x^y^ , contrary to other definitions where 
a bicubic sphne only has terms x^y^ with < i + j < 3. 
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A. 2 Spline fitting 



For any value of the parameters fk^i the coefficients Cfj - i.e. the whole spline surface S{x, y) - 
are known unambiguously. This way we consistently parameterized the spline surface (for the 
purpose of the QCD equation of state, x plays the role of the temperature and y the role of the 
quark mass ratio R). Now we want to set the spline parameters fk,i such, that derivatives of 
the spline surface are as close to some previously measured values as possible. Note that this 
way the function S{x, y) will be undetermined up to an overall constant, since the translation 
S{x, y) —7- S{x, y) + A does not influence the derivatives. This symmetry will be taken into 
account in the following. 

Let us consider N number of points and say that at each point {q^^q^) we measured 
the derivative in the x and y directions: and with errors AD^) and AD^) [m = 

. . . N — 1). Being "close" can be quantifled by minimizing 



Af-l 



m=0 




(A.5) 



Since S and thus dS/dx and dS/dy are linear in f^^i, this function has a quadratic depen- 
dence on the parameters fk,i- This enables us to search for the minimum of x^{fk,i) 



by solving a system of linear equations 



(A.6) 



M^tnjk,i = Ki,n2, ni = . . . K - 1, n2 = . . . L - 1 (A.7) 

Due to the above mentioned translational invariance of the solution, this system of equations 
is underdetermined and thus the inverse of M does not exist. This can also be seen by checking 
that M has a zero eigenvalue, corresponding to the eigenvector (1, 1,1,.. .), or, in other words, 
each row of M adds up to zero. Physically this means that one can set e.g. the flrst element 
of / to zero, i.e. leave the flrst column of M. To obtain an invertible matrix one now has 
to drop one of its rows, for example the flrst^. This way one arrives at a matrix M' of size 
{KL — 1) X [KL — 1). In the same manner we deflne V to be the vector composed from the 
last KL — 1 elements of V. One can then complement the solution {M')~^V' with a zero in 
the flrst element to obtain the flnal result which satisfles /o,o = 0. 



^It can be checked that after ehminating the first column, any row of the matrix M can be reproduced by 
a linear combination of the other rows, i.e. it is indeed correct to drop an arbitrary row. 
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In order to explicitly write the elements of the matrix M and the vector V let us analyze 
the dependence of oii the parameters fk^i- To this end we define k{m) and /(m) as the 
indices of the grid square that contains the mth measurement, i.e. 



(A.8) 



and define and rjm as the value of the dimensionless coordinate on that grid square corre- 
sponding to qm^ and (just as in (A.2)): 



{x) _ 
q-m Xk{m) 



W 



(x) 



Vn 



(y) 

Qm - yi(m) 



W 



(y) 



(A.9) 



Also, in order to be able to express dS/dx and dS/dy let us define the following matrices: 

3 



J^nun2 _ \^ (Y-l\n2 ( Y'^V^ ■ y^'V?//""^ ■ n 

ij=0 



(A.IO) 



prii,n2 
m 



i,j=0 



With the matrices E and F the expression for x in (A. 5) can be rewritten as 



N-l 



X 



m=0 



-' ^m Jni,n2 ^ 



(x)' 



m Jni,n2 ^ 



(All) 



which implies that in the system of linear equations (A. 7) to be solved appear 



Af-l 



M: 



ni,n2 



n,i,n2 



m=0 
N~l 

/ J \ m ) m m ' \ m ) in 



(A.12) 



(y) 

m 



m=0 



Using the actual form of M and V the system of linear equations in (A. 7) can be solved^ for 
fk^i- With the spline parameters the spline coefficients are also determined through (A. 4). 

Since the number of measurements is 2N and we have K ■ L — 1 independent spline parame- 
ters, the degrees of freedom of this fit is given by the difference dof = 2N — K ■ L+1. This way 
the freedom corresponding to the translational symmetry is transformed out: since /o,o = 0, 
the result S{x, y) is now set such that S{xo, yo) = holds. Note that the sphne function can 



^The system of linear equations can be solved using e.g. the Lapack library. 
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also be transformed easily to satisfy some a priori known reference equation S{xr.,yr) = Sr by 
a simple translation S — j- S + {Sr — S{xr,yr))- 

Note that the iii (A. 5) represents a situation where the input data for the derivatives 
D(^) and D'^y^ are uncorrelated. However, this is usually not the case and one has to take 
into account the correlation between the measurements and therefore the function will 
contain additional terms. Usually measurements at different m values are independent, but 
the correlation of Dx{m) and Dy{m) at the same m is not negligible, since e.g. in a lattice 
calculation these are determined using the same configurations. This leads to the following 
y2 : 

-^corr 



M~l 

\' =T 

A. corr / j 
m=0 



where Q(m) is the 2x2 correlation matrix consisting of the correlators of the two measured 
derivatives at the mth point. The corresponding system of linear equations has the same form 
as (A. 7); only now on the left and right hand side enter 



]\/Tk,l \ ^ f/O^l Tpk,l Tpni,n2 _i_ 0/0 — 1 { fpk,l rpni,n2 _i_ rpni,n2 rpk,l\ 

m=0 

m=0 

A. 3 Stable solutions 

The above prescribed method is bound to give the function S{x, y) for which the sum of 
deviations is the smallest. The solution on the other hand also depends on the number and 
the position of the nodepoints, and these have to be tuned appropriately in order to determine 
the surface that fits best. 

If the number of nodepoints K ■ L is small compared to the number of measurements N, 
then the reduced sum of deviations will be large (x^/dof ^ 1) and the spline function S{x, y) 
may not be a good approximation to the surface sought for. On the other hand if K-L is large^. 



^Obviously the inequality K ■ L — 1 < 2N should hold otherwise the problem is underdetermined. 
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the best fit will become an oscillatory function which has the correct derivatives everywhere 
(i.e. x^/dof ~ 1), but is probably not the "right" solution, especially when one knows that S 
should be monotonic^ in e.g. x. Note that this unwanted feature is a characteristic of spline 
functions even in one dimension. 




Figure A. 2: One-dimensional slice of a two-dimensional fit. Shown are the input data D^^^ and the 
derivative dS/dx of the spline surface obtained using two nodepoint sets. The two sets differ only in 
one gridpoint, see (A. 15). On the left side the number of nodepoints in the x direction is K = 10; a 
slight change in the nodepoints results in no visible change in the solution (orange and blue curves are 
on top of each other). On the other hand, on the right side K = 14, and a similar change dramatically 
distorts the solution. The corresponding values for D are 0.004 (left) and ~ 10^ (right). 



Thus we need a measure of how "right" the solution is. A useful way to define this property 
is to investigate how much / changes as the nodepoints are modified, since the oscillatory 
solutions are unstable even under a small change in the nodepoints. This way we can filter out 
the stable, realistic solutions (see illustration on figure A. 2). 

We define modified nodepoint sets x*-"^ as 

=Xk + e, if fc = a 
(q) (A.15) 
xj^ = Xk, otherwise 

with € some small number, e.g. e = {xk~i — Xq)/K/10, and in the same manner for y^^\ 
Now we carry out the fitting procedure for each of the modified nodepoint sets, resulting in 



^For example this is the case for the pressure log Z. as a function of the temperature. 
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the modified solutions f^°'^ and f^^\ The sum of the relative differences between the original 
solution / and the modified solutions 

f (/3) _ f 
Jk,l Jk,l 



1 1 



0=0 kl 



f(a) _ f 
Jkl Jk,l 



f> 



kl 



1 ^"^ 1 

/3=0 k,l 



f> 



kl 



(A.16) 



will indeed serve as an indicator of the stability of the fit. If this relative change D is under a 
few percents then the fit can be considered stable. 



A. 4 Systematics of the method 

The statistical error astat of the result from the spline fitting method described above can be 
determined using the standard jackknife algorithm, i.e. the system of linear equations (A. 7) 
needs to be solved for each jackknife sample^. The systematic error on the other hand can be 
determined by varying the nodepoints {x^, yi}. Based on experience the number of nodepoints 
may range from M/2 to M; usually an equidistant nodepoint-set can already produce a small 
value for x^/dof, but increasing the density of nodepoints in areas where the function changes 
rapidly (i.e. where the measured derivatives are large) can further help to improve the fit 
quality. 

Accordingly, a straightforward way to determine the systematic error is to generate various 
nodepoint sets with different number (and position) of gridpoints. Then, for each set r of 
the nodepoints the fit is carried out resulting in a spline function S^- and an indicator Gr = 
(Xcorr/dof) ^ of the fit quality. Results which are in the above detailed sense not stable should 
be filtered out at this point. Then at each point the systematic error a sys is determined by 

(7sys{.x,y) = sj {Srix,yy)Q - {Srix,y))l (A.17) 

with 

{Or)a = Y.'^rGr/j2Gr (A.18) 

T r 

Thus the total error can be estimated to be 

(^tot = \l(Tlys+^ltat (A. 19) 

The above discussed method can be compared to conventional integration in two dimen- 
sions [13]. Based on this comparison we conclude that the multidimensional spline algorithm 

^Note that this can be performed in a single Lapack call that solves the equations for different vectors on 
the right hand side. 
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effectively processes the input data and results in much smaller statistical errors than the usual 
integral method. Furthermore, the main advantage of the spline fitting procedure is that it 
better estimates the systematics of the result. While for ordinary integration only a single path 
in the space of the parameters is considered, the spline method is in some sense equivalent to 
taking into account all possible integration paths at the same time. Our results indicate that 
the contribution coming from these generalized paths cannot be neglected in order to estimate 
the systematics. 
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QCD thermodynamics with dynamical fermions 

Gergely Endrodi 

- Summary - 

In this thesis I studied the finite temperature transition between confined and deconfined matter 
at zero and nonzero quark densities. The findings are relevant for the understanding of the 
evolution of the early Universe and contemporary and upcoming heavy ion experiments. Our 
results were obtained using large-scale simulations using lattices with various physical extents 
and lattice spacings, at physical values of the parameters of the theory. In all three projects we 
thoroughly checked finite size effects and discretization effects, thereby controlling the two most 
important possible sources of systematic error. In order to ensure that the extrapolations lead 
to the correct continuum limit, all the necessary additive and multiplicative renormalizations 
were applied to the studied observables. 

In chapter 4 the study of the phase diagram of QCD at small chemical potentials was 
presented. We determined the curvature of the transition line which separates the chirally 
symmetric and asymmetric phases using two observables that are sensitive to this transition. 
Our results also contain information about the relative change in the width of the transition. 
This suggests that the strength of the transition decreases very slightly, therefore the weak 
crossover that is present at zero density persists up to quark chemical potentials fiu,d ^ Tc. 

In chapter 5 I showed results regarding the equation of state in full QCD. For this project 
I developed a new method to determine the pressure from the raw lattice data. This approach 
improves on the conventionally applied integral method such that it better estimates the sys- 
tematic error of the result, and is also capable of quantifying the quark mass dependence of 
the pressure. From the pressure we derived various other thermodynamic observables, which 
can be directly used to set the equation of state in e.g. hydrodynamic models. We estimate 
the effect of a charm quark in the pressure and show that it might be important already at 
T ~ 2Tc. Furthermore, we observe that the pressure, the energy density and the entropy 
density are roughly 15-20% under their Stefan-Boltzmann limits at our highest temperature 

A these values of the temperature a reliable comparison with perturbation theory cannot 
be carried out since such perturbative expansions are only convergent at much larger temper- 
atures. In chapter 6 I report results of the equation of state in pure gauge theory at extremely 
high temperatures ~ lOOOTc. We perform simulations using lattices of size up to Ng/Nt = 24 
and gather statistics that exceed results in the literature by at least an order of magnitude to 
determine the high temperature signal in the interaction measure. At large T we successfully 
match our results with perturbation theory and we are able to fit the unknown 0{g^) coeffi- 
cient of the expansion. We also quantify the non-perturbative contribution to the interaction 
measure, which is proportional to T^. 



QCD termodinamika dinamikus fermionokkal 

Endrodi Gergely 

- Osszefoglalo - 

Disszertaciomban az erosen kolcsonhato anyag bezart, illetve plazma allapota kozotti veges 
homersekleti atmenetet vizsgaltam nulla es veges kvark siirusegeknel. Az eredmenyek mind 
a korai Univerzum fejlodesenek, mind pedig napjaink nehezion-iitkoztetesi kiserleteinek szem- 
pontjabol relevansak. Az eredmenyeinket nagy szamitasigenyii racs szimulaciok hasznalataval 
kaptuk; a szimulaciokat szamos kiilonbozo meretii es racsallandoju racson vegeztiik el, az 
elmelet parametereinek fizikai erteke mellett. Mindharom bemutatott projektben megmertiik 
a veges meret effektusokat es a diszkretizacios hibakat, ezaltal iigyelve a ket legfontosabb 
lehetseges szisztematikus hibaforrasra. A vizsgalt mennyisegek minden sziikseges additiv es 
multiplikatiV renormalasat elvegeztiik. 

A 4. fejezetben a QCD fazisdiagramjanak vizsgalatat mutattam be. Meghataroztuk 
az atmeneti vonal gorbiiletet ket olyan mennyiseg segi'tsegevel, amelyek erzekenyek erre az 
atmenetre. Az eredmenyek az atmenet erossegenek relatiV valtozasarol is szamot adnak. 
Megallapitottuk, liogy az atmenet csekely mertekben gyengiil, igy a /i = 0-nal megvalosulo 
crossover alacsony kemiai potencialnal fiu,d ^ is fennmarad. 

Az 5. fejezetben a QCD allapotegyenletet vizsgaltam dinamikus fermionok jelenleteben. 
Ehhez a projekthez kifejlesztettem egy uj modszert a nyomas nyers racs adatokbol torteno 
meghatarozasara. Ez a technika a hagyomanyosan alkalmazott integral modszernel effektiveb- 
ben becsli meg a szisztematikus hibat, es a nyomas kvarktomeg-fiiggeset is kozvetleniil megadja. 
A nyomasbol szamos mas termodinamikai mennyiseget is meghataroztunk, amelyek kozvetleniil 
felhasznalliatoak az allapotegyenlet beallftasara peldaul hidrodinamikai modellekben. Megmer- 
tiik a charm kvark jarulekat a nyomashoz, es azt talaltuk, hogy mar ~ 2Tc-n jelentos jarulekrol 
van szo. Megallapitottuk tovabba, hogy a nyomas, az energiasiiriiseg es az entropiasiiriiseg 
koriilbeliil 15-20%-kal alulmiiljak Stefan-Boltzmann limeszeiket a legnagyobb homersekletiink- 
nel (T ~ 6Tc) is. 

Ezeken a homersekleteken nem lehetseges osszehasonh'tani a racs eredmenyeket a per- 
turbacioszamitassal, ugyanis a perturbatiV sorfejtesek csak rendkiViil magas homersekleten 
konvergalnak. A 6. fejezetben ezert extrem magas homersekletekig (T ~ lOOOTc) megmertiik a 
nyomast tiszta mertekelmeletben; ehhez az alkalmazott racsok teriranyii kiterjedeset rendkfviil 
nagynak (akar Ng/Nf = 24- nek) kellett megvalasztani. Igen nagy szamban generaltunk kon- 
figuraciokat; a kapott statisztikus hiba az irodalomban fellelheto eredmenyeket legalabb egy 
nagysagrenddel feliilmiilja. Magas homersekleten osszehasonlitottuk az eredmenyeket a javitott, 
illetve felosszegzett perturbacioszamftas kepleteivel, es illesztettiik az 0{g^) rend szabad egyiitt- 
hatojat a racs adatok segi'tsegevel. Vegiil ket egyszerii fiiggveny formajaban meghataroztuk a 
kolcsonhatasi mertekhez jarulo nemperturbatfv tagot, amely T^-tel aranyos. 



